if n_q==n_y			# Output weight
    Q_x =  C'*Q*C;		# Weighting on x
  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")
    A_u = A-B*k;
  else

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);

  ## Open the file 
  filename = sprintf("%s.h", name);
  fid = fopen(filename,"w");

  ## 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
      if j<N_U
	comma = ",";

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
  ## A_u         composite system matrix for U* generation 
  ##             one square matrix (A_ui) row for each system input
  ##             each A_ui generates U*' for ith system input.
  ## 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
    n_x = 0;
    n_y = 0;
    n_u = 0;
  endif
  
  no_system = n_x==0;

  [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
  
  if n_u>0
    if n_u!=length(U)/n_U

  endif
  

  [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
      j = (i-1)*n_U;
      range = j+1:j+n_U;
      Utilde = [Utilde; U(range,1)'];
    endfor

    ## Composite A matrix
    if no_system
      AA = A_u;
    else
      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;
      k = (n_u-i)*n_U;
      range = j+1:j+n_U;