Index: mttroot/mtt/lib/control/PPP/ppp_qp.m
==================================================================
--- mttroot/mtt/lib/control/PPP/ppp_qp.m
+++ mttroot/mtt/lib/control/PPP/ppp_qp.m
@@ -1,14 +1,15 @@
-function [u,U,J] = ppp_qp (x,W,J_uu,J_ux,J_uw,Us0,Gamma,gamma)
+function [u,U,iterations] = ppp_qp (x,W,J_uu,J_ux,J_uw,Us0,Gamma,gamma,mu)
 
   ## usage:  [u,U] = ppp_qp (x,W,J_uu,J_ux,J_uw,Gamma,gamma)
   ## INPUTS:
   ##      x: system state    
   ##      W: Setpoint vector
   ##      J_uu,J_ux,J_uw: Cost derivatives (see ppp_lin)
   ##      Us0: value of U* at tau=0 (see ppp_lin)
   ##      Gamma, gamma: U constrained by Gamma*U <= gamma 
+  ##      mu  Parameter of qp_mu
   ## Outputs:
   ##      u: control signal
   ##      U: control weight vector
   ##
   ## Predictive pole-placement of linear systems using quadratic programming
@@ -16,10 +17,15 @@
   ## Use ppp_lin to generate J_uu,J_ux,J_uw
   ## Use ppp_cost to evaluate resultant cost function
 
   ## Copyright (C) 1999 by Peter J. Gawthrop
   ## 	$Id$	
+
+  if nargin<9
+    mu = 0
+  endif
+
 
   ## Check the sizes
   n_x = length(x);
 
   [n_U,m_U] = size(J_uu);
@@ -31,13 +37,15 @@
   if (n != n_U)||(m != n_x)
     error("J_ux should be %ix%i not %ix%i",n_U,n_x,n,m);
   endif
 
 
-  if length(gamma)>0		# Constraints exist: do the QP algorithm 
-    U = qp(J_uu,(J_ux*x - J_uw*W),Gamma,gamma); # QP solution for weights U
-    #U = pd_lcp04(J_uu,(J_ux*x - J_uw*W),Gamma,gamma); # QP solution for weights U
+  if length(gamma)>0		# Constraints exist: do the QP algorithm
+    [U,iterations] = qp_mu(J_uu,(J_ux*x - J_uw*W),Gamma,gamma,mu); # QP solution for weights U	
+
+    ##U = qp(J_uu,(J_ux*x - J_uw*W),Gamma,gamma); # QP solution for weights U
+    ##U = pd_lcp04(J_uu,(J_ux*x - J_uw*W),Gamma,gamma); # QP solution for weights U
     u = Us0*U;			# Control signal
   else			# Do the unconstrained solution
     ## Compute the open-loop gains
     K_w = J_uu\J_uw;
     K_x = J_uu\J_ux;

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,6 @@
-function [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,Delta_ol,movie)
+function [T,y,u,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
@@ -12,12 +12,17 @@
   if nargin<19			# No intermittent control
     Delta_ol = 0;
   endif
 
   if nargin<20			# No movie
+    mu = 0;
+  endif
+
+  if nargin<21			# No movie
     movie = 0;
   endif
+
 
   ## Check some sizes
   [n_x,n_u,n_y] = abcddim(A,B,C,D);
 
   [n_x0,m_x0] = size(x_0);
@@ -47,11 +52,11 @@
 
   ## Set up various time vectors
   dt = t(2)-t(1);		# Time increment
 
   ## Make sure Delta_ol is multiple of dt
-  Delta_ol = floor(Delta_ol/dt)*dt
+  Delta_ol = floor(Delta_ol/dt)*dt;
 
   if Delta_ol>0			# Intermittent control
     T_ol = 0:dt:Delta_ol-dt;	# Create the open-loop time vector
   else
     T_ol = 0;
@@ -73,10 +78,11 @@
   x = x_0;			# Initialise state
 
   ## Initialise the saved variable arrays
   X = [];
   u = [];
+  Iterations = [];
   du = [];
   J = [];
   tick= time;
   i = 0;
   disp("Simulating ...");
@@ -88,11 +94,11 @@
     ## Composite constraints 
     Gamma = [Gamma_u; Gamma_y];
     gamma = [gamma_u; gamma_y];
     
     ## Compute U(t)
-    [uu U] = ppp_qp (x,W,J_uu,J_ux,J_uw,Us0,Gamma,gamma); # 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];
 
@@ -107,10 +113,11 @@
 
       ## Simulate the system
       i = i+1;
       X = [X x];		# Save state
       u = [u ut];		# Save input
+      Iterations = [Iterations iterations]; # Save iteration count
       x = Ad*x + Bd*ut;	# System
 
 #       if movie			# Plot the moving horizon
 # 	tau = T(1:n_T-i);	# Tau with moving horizon
 # 	tauT = T(i+1:n_T);	# Tau with moving horizon + real time
@@ -121,17 +128,17 @@
   endfor
   
   ## Save the last values
   X = [X x];		# Save state
   u = [u ut];		# Save input
+  Iterations = [Iterations iterations]; # Save iteration count
 
   tock = time;
-  Iterations = length(T_cl)
   Elapsed_Time = tock-tick
   y = C*X + D*u;		# System output
 
   T = 0:dt:t+Delta_ol;		# Overall time vector
 
 endfunction