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[,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,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
error("U must be a column vector with n_u*n_U components");
endif
else
n_u = length(U)/n_U; # Deduce n_u from U if no system
endif
[n_x0,m_x0] = size(x_0);
if n_x0!=n_x
error(sprintf("x_0 must be a column with length %i", n_x));
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;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;
Utilde = [Utilde; zeros(1,j), U(range,1)', zeros(1,k)];
endfor
## Create the full A_u matrix (AA_u) with the A_i s on the diagonal
# AA_u = [];
# for i = 1:n_u
# AA_u = ppp_aug(AA_u,ppp_extract(A_u,i));
# endfor
AA_u = ppp_inflate(A_u);
## Composite A matrix
if no_system
AA = AA_u;
else
Z = zeros(n_U*n_u,n_x);
AA = [A B*Utilde
Z AA_u];
endif
xx_0 = [x_0;ones(n_U*n_u,1)]; # Composite initial condition
endif
## Initialise
xs = []; # x star
xu = []; # x star
ys = []; # y star
us = []; # u star
n_xx = length(xx_0); # Length of composite state
## Compute the star variables
for t=tau
xxt = expm(AA*t)*xx_0; # Composite state
xst = xxt(1:n_x); # x star
xut = xxt(n_x+1:n_xx); # x star
if length(C)>0
yst = C*xst; # y star
else
yst = [];
endif
ust = Utilde*xut; # u star
xs = [xs xst]; # x star
xu = [xu xut]; # x star
ys = [ys yst]; # y star
us = [us ust]; # u star
endfor
endfunction