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
@@ -19,16 +19,39 @@
   if nargin<1			# Default name to dir name
     names = split(pwd,"/");
     [n_name,m_name] = size(names);
     Name = deblank(names(n_name,:));
   endif
+
+  if nargin<6
+    p_c.N = 50;
+  endif
+
+  if nargin<7
+    p_o.sigma = 1e-1;
+  endif
 
   ## System
   sys = mtt2sys(Name);		# Create system
-  [A,B,C,D] = sys2ss(sys);	# SS form
+  [A,B,C_0,D_0] = sys2ss(sys);	# SS form
+
+  ## Extract matrices for controlled and constrained outputs.
+  if !struct_contains(p_c,"I_0") # Indices for controlled outputs
+    p_c.I_0 = 1:n_y
+  endif
+  if !struct_contains(p_c,"I_1") # Indices for constarined outputs
+    p_c.I_1 = 1:n_y
+  endif
+
+  C = C_0(p_c.I_0,:)
+  C_c = C_0(p_c.I_1,:);
+  D = D_0(p_c.I_0,:)
+  D_c = D_0(p_c.I_1,:);
   [n_x, n_u, n_y] = abcddim(A,B,C,D); # Dimensions
-  
+  [n_x, n_u, n_y_c] = abcddim(A,B,C_c,D_c); # Dimensions
+
+
   if nargin<2
     Simulate = 1;
   endif
   
   if nargin<3
@@ -41,17 +64,10 @@
   
   if nargin<5
     x_0 = zeros(n_x,1);
   endif
   
-  if nargin<6
-    p_c.N = 50;
-  endif
-
-  if nargin<7
-    p_o.sigma = 1e-1;
-  endif
 
   if !struct_contains(p_c,"delta_ol")
     p_c.delta_ol = 0.5;	# OL sample interval
   endif
   
@@ -134,11 +150,10 @@
     ##p_o.method = "continuous";
     p_o.method = "intermittent";
   endif
   
 
-
   ## Check w.
   [n_w,m_w] = size(w);
   if ( (n_w!=n_y) || (m_w!=1) )
     error(sprintf("ppp_lin_run: w must a column vector with %i elements",n_y));
   endif
@@ -174,19 +189,17 @@
     if (p_c.n_U!=M_u)
       error("A_u must be square");
     endif
 
     ## Checks
-    [ol_zeros, ol_poles] = sys2zp(sys)
     cl_poles = eig(A - B*k_x)
     t_max = 1/min(abs(cl_poles))
     t_min = 1/max(abs(cl_poles))
   endif
-
+  
   ## Initial control U
   U = zeros(p_c.n_U,1)	
-
 
   ## Short sample interval
   dt = p_c.delta_ol/p_c.N;
 
   ## Observer design
@@ -279,11 +292,11 @@
 	[Gamma_u, gamma_u] = \
 	    ppp_input_constraints(p_c.A_u,p_c.Tau_u,p_c.Min_u,p_c.Max_u);
 	
 	## Output constraints
 	[Gamma_y,gamma_y] = \
-	    ppp_output_constraints(A,B,C,D,x_est,p_c.A_u,\
+	    ppp_output_constraints(A,B,C_c,D_c,x_est,p_c.A_u,\
 				   p_c.Tau_y,p_c.Min_y,p_c.Max_y);
 	
 	## Composite constraints - t=0
 	Gamma = [Gamma_u; Gamma_y];
 	gamma = [gamma_u; gamma_y];