function [t,y,u,X_est,y_c,t_e,y_e,e_e,p_c,p_o] = ppp_lin_run (Name,Simulate,ControlType,w,x_0,p_c,p_o)
## usage: [t,y,u,y_c,t_e,y_e,e_e,p_c,p_o] = ppp_lin_run (Name,Simulate,ControlType,w,x_0,p_c,p_o)
##
##
## Linear closed-loop PPP of lego system (and simulation)
##
## Name: Name of system (in mtt terms)
## Simulate = 0: real thing
## Simulate = 1: simulate
## Control = 0: step test
## Control = 1: PPP open-loop
## Control = 2: PPP closed-loop
## w is the (constant) setpoint
## par_control and par_observer are structures containing parameters
## for the observer and controller
##Defaults
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_0,D_0] = sys2ss(sys); # SS form
[n_x, n_u, n_y] = abcddim(A,B,C_0,D_0);
## 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 constrained 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
ControlType = 2;
endif
if nargin<4
w = ones(n_y,1);;
endif
if nargin<5
x_0 = zeros(n_x,1);
endif
if !struct_contains(p_c,"delta_ol")
p_c.delta_ol = 0.5; # OL sample interval
endif
if !struct_contains(p_c,"T")
p_c.T = 10; # Last time point.
endif
if !struct_contains(p_c,"Iterations")
p_c.Iterations = 5; # Number of interations, total =T*Iterations
endif
if !struct_contains(p_c,"augment")
p_c.augment = 0; # Augment basis funs with constant
endif
if !struct_contains(p_c,"integrate")
p_c.integrate = 0;
endif
if !struct_contains(p_c,"Tau_u")
p_c.Tau_u = [];
p_c.Min_u = [];
p_c.Max_u = [];
endif
if !struct_contains(p_c,"Tau_y")
p_c.Tau_y = [];
p_c.Min_y = [];
p_c.Max_y = [];
endif
if !struct_contains(p_c,"Method")
p_c.Method = "lq";
endif
if struct_contains(p_c,"Method")
if strcmp(p_c.Method,"lq")
p_c.Q = eye(n_y);
p_c.n_U = n_x;
if !struct_contains(p_c,"R")
p_c.R = (0.1^2)*eye(n_u);
endif
elseif strcmp(p_c.Method,"original");
if !struct_contains(p_c,"A_w")
p_c.A_w = 0;
endif
if !struct_contains(p_c,"A_u")
p_c.n_U = n_x;
a_u = 2.0;
p_c.A_u = laguerre_matrix(p_c.n_U,a_u);
if p_c.augment # Put in constant term
p_c.A_u = ppp_aug(0,p_c.A_u);
endif
endif
else
error(sprintf("Method %s not recognised", p_c.Method));
endif
endif
if !struct_contains(p_c,"tau") # Time horizon
if strcmp(p_c.Method,"lq")
p_c.tau = [0:0.1:1]*2;
elseif strcmp(p_c.Method,"original");
p_c.tau = [10:0.1:11];
else
error(sprintf("Method %s not recognised", p_c.Method));
endif
endif
if !struct_contains(p_c,"A_e")
p_c.A_e = []; # No extra modes
endif
if !struct_contains(p_o,"x_0")
p_o.x_0 = zeros(n_x,1);
endif
if !struct_contains(p_o,"method")
##p_o.method = "continuous";
##p_o.method = "intermittent";
p_o.method = "remote";
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
## Initialise
x_est = p_o.x_0;
## Initialise simulation state
x = x_0;
y_i = C*x_0
if ControlType==0 # Step input
I = 1; # 1 large sample
p_c.delta_ol = p_c.T; # I
K_w = zeros(p_c.n_U,n_y);
K_w(1,1) = 1;
K_w(2,1) = -1;
K_x = zeros(p_c.n_U,n_x);
else
I = ceil(p_c.T/p_c.delta_ol) # Number of large samples
if strcmp(p_c.Method, "original")
[k_x,k_w,K_x,K_w,Us0,J_uu,J_ux,J_uw,J_xx,J_xw,J_ww] =\
ppp_lin(A,B,C,D,p_c.A_u,p_c.A_w,p_c.tau); # Design
elseif strcmp(p_c.Method, "lq") # LQ design
[k_x,k_w,K_x,K_w,Us0,J_uu,J_ux,J_uw,J_xx,J_xw,J_ww,A_u] \
= ppp_lin_quad (A,B,C,D,p_c.tau,p_c.Q,p_c.R,p_c.A_e);
p_c.A_u = A_u;
else
error(sprintf("Control method %s not recognised", p_c.Method));
endif
##Sanity check A_u
[p_c.n_U,M_u] = size(p_c.A_u);
if (p_c.n_U!=M_u)
error("A_u must be square");
endif
## Checks
cl_poles = eig(A - B*k_x)
ol_poles = eig(A)
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
G = eye(n_x); # State noise gain
sigma_x = eye(n_x); # State noise variance
Sigma = p_o.sigma*eye(n_y) # Measurement noise variance
if strcmp(p_o.method, "intermittent")
Ad = expm(A*p_c.delta_ol); # Discrete-time transition matrix
if (ControlType==2) #
[L, M, P, obs_poles] = dlqe(Ad,G,C,sigma_x,Sigma);
else
L = zeros(n_x,n_y);
obs_poles = eig(Ad);
endif
elseif strcmp(p_o.method, "continuous")
Ad = expm(A*dt); # Discrete-time transition matrix
A_ud = expm(p_c.A_u*dt); # Discrete-time input transition
if (ControlType==2) #
[L, M, P, obs_poles] = dlqe(Ad,G,C,sigma_x,Sigma);
else
L = zeros(n_x,n_y);
obs_poles = eig(Ad);
endif
elseif strcmp(p_o.method, "remote")
L = zeros(n_x,n_y);
obs_poles = [];
else
error(sprintf("Observer method ""%s"" unknown", p_o.method));
endif
## Display the poles
obs_poles
## Write the include file for the real-time function
## Use double length to allow for overuns
overrun = 2;
Ustar = ppp_ustar (p_c.A_u, n_u, [0:dt:overrun*p_c.delta_ol], 0,0);
if p_c.integrate # Integrate Ustar
disp("Integrating Ustar");
Ustar = cumsum(Ustar)*dt;
endif
disp("Writing Ustar.h ...");
ppp_ustar2h(Ustar);
disp("done.");
## Control loop
y = [];
y_c = [];
u = [];
t = [];
y_e = [];
X_est = [];
t_e = [];
e_e = [];
tick = time;
i=0;
for j=1:p_c.Iterations
for k=1:I
tim=time; # Timing
i++;
if Simulate # Exact simulation
X = x; # Current (simulated) state
else # The real thing
if strcmp(p_o.method, "remote")
[t_i,y_i,u_i,X] = ppp_put_get_X(U); # Remote-state interface
else
[t_i,y_i,u_i] = ppp_put_get(U); # Generic interface to real-time
endif
endif
## Observer
if strcmp(p_o.method, "intermittent")
[x_est y_est y_new, e_est] = ppp_int_obs \
(x_est,y_i,U,A,B,C,D,p_c.A_u,p_c.delta_ol,L);
elseif strcmp(p_o.method, "continuous")
Ui = U; # U at sub intervals
for k = 1:p_c.N
[x_est y_est y_new e_est] = ppp_int_obs \
(x_est,yi(:,k),Ui,A,B,C,D,p_c.A_u,dt,L);
Ui = A_ud'*Ui;
y_e = [y_e; y_new'];
e_e = [e_e; e_est'];
endfor
elseif strcmp(p_o.method, "remote")
## predict from remote state (with zero L)
if (ControlType==2) # Closed-loop
[x_est y_est y_new e_est] = ppp_int_obs \
(X,y_i,U,A,B,C,D,p_c.A_u,p_c.delta_ol,zeros(n_x,1));
## x_est = X; y_est=y_i; y_new=y_i; e_est=0;
else # Open-loop
[x_est y_est y_new e_est] = ppp_int_obs \
(x_est,y_i,U,A,B,C,D,p_c.A_u,p_c.delta_ol,zeros(n_x,1));
endif
endif
## Simulation (based on U_i)
if Simulate
t_sim = [1:p_c.N]*dt; # Simulation time points
[yi,ui,xsi] = ppp_ystar(A,B,C,D,x,p_c.A_u,U,t_sim); # Simulate
x = xsi(:,p_c.N); # NEXT state
ti = [(i-1)*p_c.N:i*p_c.N-1]*dt;
y_i = yi(1); # Current output
t_i = ti(1);
endif
##Control
if ( length(p_c.Tau_u)==0&&length(p_c.Tau_y)==0 )
U = K_w*w - K_x*x_est;
else
## Input constraints
[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_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];
[u_qp,U,n_active] = ppp_qp \
(x_est,w,J_uu,J_ux,J_uw,Us0,Gamma,gamma,1e-6,1);
endif
## Allow for the delay
##U = expm(p_c.delta_ol*p_c.A_u)*U;
## Save data
if Simulate
t = [t;ti'];
y = [y;yi'];
X_est = [X_est;x_est'];
y_c = [y_c;(C_c*xsi)'];
u = [u;ui'];
else
t = [t;t_i];
y = [y;y_i'];
X_est = [X_est;x_est'];
u = [u;u_i'];
endif
if strcmp(p_o.method, "intermittent")||strcmp(p_o.method, "remote")
y_e = [y_e; y_new'];
e_e = [e_e; e_est'];
t_e = [t_e; t_i];
endif
if !Simulate
delta_comp = time-tim;
usleep(floor(1e6*(p_c.delta_ol-delta_comp-0.01)));
endif
endfor # Main loop
w = -w;
endfor # Outer loop
if !Simulate
ppp_put_get(0*U); # Reset to zero
endif
if strcmp(p_o.method, "continuous")
t_e = t;
endif
average_ol_sample_interval = (time-tick)/i
## Put data on file (so can use for identification)
filename = sprintf("%s_ident_data.dat",Name);
eval(sprintf("save -ascii %s t y u",filename));
endfunction