Index: mttroot/mtt/lib/control/PPP/ppp_are.m
==================================================================
--- mttroot/mtt/lib/control/PPP/ppp_are.m
+++ mttroot/mtt/lib/control/PPP/ppp_are.m
@@ -26,12 +26,12 @@
   elseif n_q==n_x		# State weight
     Q_x = Q;
   else
     error(sprintf("Q (%ix%i) must be %ix%i or %ix%i",n_q,n_q,n_y,n_y,n_x,n_x));
   endif
-  Q_x
-  [k, P, poles] = lqr (A, B, Q_x, R) # Algebraic Riccati solution
+
+  [k, P, poles] = lqr (A, B, Q_x, R); # Algebraic Riccati solution
 
   ## Basis functions
   if strcmp(A_type,"companion")
     A_u = compan(poly(poles));
   elseif strcmp(A_type,"feedback")

Index: mttroot/mtt/lib/control/PPP/ppp_ustar2h.m
==================================================================
--- mttroot/mtt/lib/control/PPP/ppp_ustar2h.m
+++ mttroot/mtt/lib/control/PPP/ppp_ustar2h.m
@@ -1,12 +1,16 @@
-function ppp_ustar2h (Ustar,name)
+function ppp_ustar2h (Ustar,DT,name)
 
   ## usage:  ppp_Ustar2h (Ustar[,name])
   ##
   ## 
 
   if nargin<2
+    DT = 1;
+  endif
+  
+  if nargin<3
     name = "Ustar";
   endif
 
   [N,N_U] = size(Ustar);
 
@@ -16,11 +20,13 @@
 
   ## Header
   header = sprintf("/*\n File %s generated by ppp_ustar2h on %s */\n", \
 		   filename, ctime(time));
 
-  def = sprintf("#define N_U %i\n#define N_T %i\n", N_U, N);
+  def = sprintf("#define N_U %i\n#define N_T %i\n#define DT %g\n", \
+		 N_U, N, DT);
+  def = sprintf("%sdouble U[N_U];\n",def);
 
   fprintf(fid, "%s%sdouble %s[N_T][N_U] = {\n",header,def,name);
   for i=1:N
     fprintf(fid, "{");
     for j=1:N_U

Index: mttroot/mtt/lib/control/PPP/ppp_ystar.m
==================================================================
--- mttroot/mtt/lib/control/PPP/ppp_ystar.m
+++ mttroot/mtt/lib/control/PPP/ppp_ystar.m
@@ -1,8 +1,8 @@
-function [ys,us,xs,xu,AA] = ppp_ystar (A,B,C,D,x_0,A_u,U,tau)
+function [ys,us,xs,xu,AA] = ppp_ystar (A,B,C,D,x_0,A_u,U,tau,Us0)
 
-  ## usage:  [ys,us,xs,xu,AA] = ppp_ystar (A,B,C,D,x_0,A_u,U,tau)
+  ## usage:  [ys,us,xs,xu,AA] = ppp_ystar (A,B,C,D,x_0,A_u,U,tau[,Us0])
   ##
   ## Computes open-loop moving horizon variables at time tau
   ## Inputs:
   ## A,B,C,D     System matrices
   ## x_0         Initial state
@@ -12,18 +12,20 @@
   ## OR
   ## A_u         square system matrix for U* generation 
   ##             same square matrix for each system input
   ## U           Column vector of optimisation coefficients  
   ## tau         Row vector of times at which outputs are computed
+  ## Us0         Initial value of U* (default ones(NU,1))
   ## Outputs:
   ## ys          y*, one column for each time tau 
   ## us          u*, one column for each time tau 
   ## xs          x*, one column for each time tau 
   ## xu          x_u, one column for each time tau 
   ## AA          The composite system matrix
   
-  ## Copyright (C) 1999 by Peter J. Gawthrop
+
+  ## Copyright (C) 1999,2005 by Peter J. Gawthrop
   ## 	$Id$	
 
   if (size(A)>0)
     [n_x,n_u,n_y] = abcddim(A,B,C,D); # System dimensions
   else
@@ -36,11 +38,10 @@
 
   [n,m] = size(A_u);		# Size of composite A_u matrix
   square = (n==m);		# Is A_u square?
   n_U = m;			# functions per input
 
-  
   [n,m] = size(U);
   if (m != 1)
     error("U must be a column vector");
   endif
   
@@ -60,10 +61,20 @@
 
   [n,m]=size(tau);
   if (n != 1 )
     error("tau must be a row vector of times");
   endif
+
+  if nargin<9
+    Us0 = ones(1,n_U);
+  endif
+  
+  [n_Us0,m_Us0] = size(Us0);
+  if (n_Us0>1)||(n_Us0>m_Us0)
+    error(sprintf("Us0 must be a row vector, not %ix%i ",n_Us0,m_Us0));
+  endif
+	  
   
   if square			# Then same A_u for each input
     ## Reorganise vector U into matrix Utilde  
     Utilde = [];
     for i=1:n_u
@@ -79,11 +90,11 @@
       Z = zeros(n_U,n_x);
       AA = [A   B*Utilde
 	    Z   A_u];
     endif
     
-    xx_0 = [x_0;ones(n_U,1)];	# Composite initial condition
+    xx_0 = [x_0;Us0'];		# Composite initial condition
   else				# Different A_u on each input
     ## Reorganise vector U into matrix Utilde  
     Utilde = [];
     for i=1:n_u
       j = (i-1)*n_U;