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