function [y,ystar,t,k_x,k_w,K_x,K_w,Us0,J_uu,J_ux,J_uw,J_xx,J_xw,J_ww,y_u,cond_uu] =  ppp_lin_sim(A,B,C,D,A_u,A_w,tau,Q,W,x_0)

  ## usage:   [ol_poles,cl_poles,ol_zeros,cl_zeros,k_x,k_w,K_x,K_w] = ppp_lin_plot (A,B,C,D,A_u,A_w,tau,Q,W,x_0)
  ##
  ## Linear PPP (Predictive pole-placement) computation with plotting
  ## INPUTS:
  ## A,B,C,D: system matrices
  ## 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.
  ## A_w: composite system matrix for W* generation 
  ##      one square matrix (A_wi) row for each system output
  ##      each A_wi generates W*' for ith system output.
  ## tau: row vector of times for optimisation (equispaced in time)
  ## Q: column vector of output weights (defaults to unity)
  ## W: Constant setpoint vector (one element per output)
  ## x_0: Initial state
  ## OUTPUTS:
  ## y : closed-loop output
  ## ystar : open-loop moving-horizon output
  ## t : time axis

  ## Copyright (C) 2001 by Peter J. Gawthrop
  ## 	$id: ppp_lin_plot.m,v 1.13 2001/01/26 16:03:13 peterg Exp $	

  ## Some dimensions
  [n_x,n_u,n_y] = abcddim(A,B,C,D);
  [n_U,m_U]=size(A_u);
  square = (n_U==m_U);		# Its a square matrix so same U* on each input
  [n_W,m_W]=size(A_w);
  if n_W==m_W			# A_w square
    n_W = n_W*n_y;		# Total W functions
  endif
  

  [n_tau,m_tau] = size(tau);

  ## Default Q
  if nargin<8
    Q = ones(n_y,1);
  endif

  ## Default W
  if nargin<9
    W = ones(n_W,1)
  endif

  ## Default x_0
  if nargin<10
    x_0 = zeros(n_x,1);
  endif

  ## Check some dimensions
  [n_Q,m_Q] = size(Q);
  if ((m_Q != 1)&&(m_Q != m_tau)) || (n_Q != n_y)
    error("Q must be a column vector with one row per system output");
  endif

  [n,m] = size(W);
  if ((m != 1) || (n != n_W))
    error("W must be a column vector with one element per system output");
  endif

  [n,m] = size(x_0);
  if ((m != 1) || (n != n_x))
    error("x_0 must be a column vector with one element per system state");
  endif

  ## Control design
  disp("Designing controller");
  [k_x,k_w,K_x,K_w,Us0,J_uu,J_ux,J_uw,J_xx,J_xw,J_ww,y_u,cond_uu] = ppp_lin  (A,B,C,D,A_u,A_w,tau,Q);

  ## Set up simulation times
  dtau = tau(2) -tau(1);		# Time increment
  t =  0:dtau:tau(length(tau));	# Time starting at zero

  ## Compute the OL step response
  disp("Computing OL response");
  U = K_w*W - K_x*x_0;
  ystar = ppp_ystar (A,B,C,D,x_0,A_u,U,t);

  ## Compute the CL step response
  disp("Computing CL response");
  y = ppp_sm2sr(A-B*k_x, B, C, D, t, k_w*W, x_0); # Compute Closed-loop control


endfunction