Index: mttroot/mtt/lib/control/PPP/ppp_qp_sim.m
==================================================================
--- mttroot/mtt/lib/control/PPP/ppp_qp_sim.m
+++ mttroot/mtt/lib/control/PPP/ppp_qp_sim.m
@@ -1,6 +1,9 @@
-function [T,y,u,X,Iterations] = ppp_qp_sim (A,B,C,D,A_u,A_w,t,Q, Tau_u,Min_u,Max_u,Order_u, Tau_y,Min_y,Max_y,Order_y, W,x_0,Delta_ol,mu,movie)
+function [T,y,u,X,Iterations] = ppp_qp_sim (A,B,C,D,A_u,A_w,t,Q,\
+					    Tau_u,Min_u,Max_u,Order_u, \
+					    Tau_y,Min_y,Max_y,Order_y, \
+					    W,x_0,Delta_ol,mu,movie)
 
   ## usage: [T,y,u,J] = ppp_qp_sim (A,B,C,D,A_u,A_w,t,Q, Tau_u,Min_u,Max_u,Order_u, Tau_y,Min_y,Max_y,Order_y, W,x_0,movie)
   ## Needs documentation - see ppp_ex11 for example of use.
   ## OUTPUTS
   ## T: Time vector
@@ -28,21 +31,15 @@
   [n_x0,m_x0] = size(x_0);
   if (n_x0 != n_x)||(m_x0 != 1)
     error(sprintf("Initial state x_0 must be %ix1 not %ix%i",n_x,n_x0,m_x0));
   endif
   
-  ## Input constraints (assume same on all inputs)
-  Gamma_u=[];
-  gamma_u=[];
-  for i=1:n_u
-    [Gamma_i,gamma_i] = ppp_input_constraint (A_u,Tau_u,Min_u,Max_u,Order_u,i,n_u);
-    Gamma_u = [Gamma_u; Gamma_i];
-    gamma_u = [gamma_u; gamma_i];
-  endfor
+  ## Input constraints 
+  [Gamma_u, gamma_u] = ppp_input_constraints(A_u,Tau_u,Min_u,Max_u);
 
   ## Output constraints
-  [Gamma_y,gamma_y] = ppp_output_constraint  (A,B,C,D,x_0,A_u,Tau_y,Min_y,Max_y,Order_y);
+  [Gamma_y,gamma_y] = ppp_output_constraints(A,B,C,D,x_0,A_u,Tau_y,Min_y,Max_y,Order_y);
 
   ## Composite constraints - t=0
   Gamma = [Gamma_u; Gamma_y];
   gamma = [gamma_u; gamma_y];
 
@@ -60,15 +57,29 @@
     T_ol = 0:dt:Delta_ol;	# Create the open-loop time vector
   else
     T_ol = [0,dt];
     Delta_ol = dt;
   endif
-  T_cl = 0:Delta_ol:t(length(t))-Delta_ol; # Closed-loop time vector
+  t_last = 2*t(length(t));
+  T_cl = 0:Delta_ol:t_last-Delta_ol; # Closed-loop time vector
+  T = 0:dt:t_last;		# Overall time vector
+ 
+  ## Lengths thereof
   n_Tcl = length(T_cl);
   n_ol = length(T_ol);
-  
+  n_T = length(T);
+
+  ## Expand W with constant last value or truncate
+  [n_W,m_W] = size(W);
+
+  if m_W>n_T
+    W = W(:,1:n_T);
+  else
+    W = [W W(:,m_W)*ones(1,n_T-m_W+1)];
+  endif
 
+  ## Compute U*
   Ustar_ol = ppp_ustar(A_u,n_u,T_ol); # U* in the open-loop interval
 
   [n,m] = size(Ustar_ol);
   n_U = m/length(T_ol);		# Determine size of each Ustar
 
@@ -84,23 +95,27 @@
   u = [];
   Iterations = [];
   du = [];
   J = [];
   tick= time;
-  i = 0;
+
   ## disp("Simulating ...");
   for t=T_cl			# Outer loop at Delta_ol
+
     ##disp(sprintf("Time %g", t));
     ## Output constraints
-    [Gamma_y,gamma_y] = ppp_output_constraint  (A,B,C,D,x,A_u,Tau_y,Min_y,Max_y,Order_y);
+    [Gamma_y,gamma_y] = ppp_output_constraints  (A,B,C,D,x,A_u,Tau_y,Min_y,Max_y,Order_y);
     
     ## Composite constraints 
     Gamma = [Gamma_u; Gamma_y];
     gamma = [gamma_u; gamma_y];
+    
+    ## Current Setpoint value
+    w = W(:,floor(t/dt)+1);
     
     ## Compute U(t) via QP optimisation
-    [uu, U, iterations] = ppp_qp (x,W,J_uu,J_ux,J_uw,Us0,Gamma,gamma,mu); # Compute U
+    [uu, U, iterations] = ppp_qp (x,w,J_uu,J_ux,J_uw,Us0,Gamma,gamma,mu); # Compute U
 
     ## Compute the cost (not necessary but maybe interesting)
 #    [J_t] = ppp_cost (U,x,W,J_uu,J_ux,J_uw,J_xx,J_xw,J_ww); # cost
 #    J = [J J_t];
 
@@ -125,11 +140,10 @@
 
   tock = time;
   Elapsed_Time = tock-tick;
   y = C*X + D*u;		# System output
 
-  T = 0:dt:t+Delta_ol;		# Overall time vector
 
 endfunction