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
@@ -1,22 +1,37 @@
-function [P,A_u,A_w,k] = ppp_are (A,B,C,D,Q,R)
+function [P,A_u,A_w,k] = ppp_are (A,B,C,D,Q,R,A_type)
 
-  ## usage:  [P,A_u,A_w] = ppp_are (A,B,C,D,Q,R)
+  ## usage:  [P,A_u,A_w] = ppp_are (A,B,C,D,Q,R,A_type)
   ##
   ## 
 
+  if nargin<1
+    disp("usage:  [P,A_u,A_w] = ppp_are (A,B,C,D,Q,R,A_type)");
+    return
+  endif
+  
+  if nargin<7
+    A_type = "feedback";
+  endif
+  
 
   ## Steady-state Linear Quadratic solution
   ## using Algebraic Riccati equation (ARE)
   Q_x =  C'*Q*C;			# Weighting on x
   [k, P, poles] = lqr (A, B, Q_x, R); # Algebraic Riccati solution
 
   ## Basis functions
-  A_u = compan(poly(poles));
-
-  ## Avoid spurious imag parts due to rounding
-  A_u = real(A_u);		
+  if strcmp(A_type,"companion")
+    A_u = compan(poly(poles));
+  elseif strcmp(A_type,"feedback")
+    A_u = A-B*k;
+  else
+    error(sprintf("A_type must be %s, not %s", "companion or feedback", A_type));
+  endif
+  
+    ## Avoid spurious imag parts due to rounding
+    A_u = real(A_u);		
 
   ## Setpoint basis functions
   A_w = 0;
 
 endfunction

Index: mttroot/mtt/lib/control/PPP/ppp_lin_run.m
==================================================================
--- mttroot/mtt/lib/control/PPP/ppp_lin_run.m
+++ mttroot/mtt/lib/control/PPP/ppp_lin_run.m
@@ -13,12 +13,10 @@
   ## Control = 2:  PPP closed-loop
   ## w is the (constant) setpoint
   ## par_control and par_observer are structures containing parameters
   ## for the observer and controller
 
-  disp("Special version");
-
   ##Defaults
   if nargin<1			# Default name to dir name
     names = split(pwd,"/");
     [n_name,m_name] = size(names);
     Name = deblank(names(n_name,:));

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,6 +1,6 @@
-function [u,U,iterations] = ppp_qp (x,W,J_uu,J_ux,J_uw,Us0,Gamma,gamma,mu,test)
+function [u,U,n_active] = ppp_qp (x,W,J_uu,J_ux,J_uw,Us0,Gamma,gamma,mu,test)
 
   ## usage:  [u,U] = ppp_qp (x,W,J_uu,J_ux,J_uw,Gamma,gamma)
   ## INPUTS:
   ##      x: system state    
   ##      W: Setpoint vector
@@ -43,22 +43,24 @@
   endif
 
 
   if length(gamma)>0		# Constraints exist: do the QP algorithm
     ## QP solution for weights U	
-    [U,iterations] = qp_mu(J_uu,(J_ux*x - J_uw*W),Gamma,gamma,mu,[],[],0,test);
+##    [U,iterations] = qp_mu(J_uu,(J_ux*x-J_uw*W),Gamma,gamma,mu,[],[],0,test);
+    [U,n_active] = qp_hild(J_uu,(J_ux*x - J_uw*W),Gamma,gamma);	# 
+##    iterations = 0;
 
     ##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
-    iterations = 0;
+    n_active = 0;
     K_w = J_uu\J_uw;
     K_x = J_uu\J_ux;
 
     ## Closed-loop control
     U = K_w*W - K_x*x;		# Basis functions weights - U(t)
     u = Us0*U;			# Control u(t)
   endif
 
 endfunction