% Adaptive MV algorithm, sec. 5.11.2.
%
% The program simulates control system with adaptive
% minimumvariance controller. The plant and the controller
% are simulated as discrete in time.
%
% The plant parameters are: polynomials A, B and delay-time
% k_delay of the transfer function z^(-k_delay) B/A and
% standard deviation sd of white noise of normal distribution.
% Parameters of RLS are: initial values of diagonal of
% covariance matrix P_init and forgetting factor alpha.
%
% The controller has the form u(i) = - (S/R) y(i) where
% S of order nS and R of order nR are polynomials in z^(-1).
% The structure of the controller is calculated according
% to the plant structure but orders nS and nR as well as
% initial values of the polynomials can be changed by the user.
% In predefined number of sample time ch_time new parameters
% of the plant can be introduced what provides time-varying
% plant. Controller parameters are identified on-line.
% Values of polynomial R and S parameters are displayed after
% simulation finishes.
%
% Number of sample times to be simulated - end_time - has also
% to be chosen by the user.
%
%---------------------------------------------------------------
% copyright 1995, by J. Moscinski and Z. Ogonowski. For use with
% the book "Advanced Control with MATLAB and SIMULINK"
% (Prentice-Hall, 1995).
%---------------------------------------------------------------
par_R = [];
par_S = [];
z = menu('Parameters of RLS', ...
'P_init = 1e5; alpha = 1', ...
'Change the parameters');
disp('')
if z == 1
P_init = 100000;
alpha = 1;
else
P_init = input('P_init = ');
alpha = input('aplha = ');
end
z = menu('Initial parameters of the plant', ...
'A = [1 -0.7555 0.0498]; B = [0.2134 0.081]; k = 1; sd = 1', ...
'Change the parameters');
disp('')
if z == 1
k_delay = 1;
A = [-0.7555 0.0498];
B = [zeros(1,k_delay-1), 0.2134 0.081];
sd = 1;
else
A =input('A = ');
A = A(2:length(A));
B = input('B = ');
B = [zeros(1,k_delay-1), B];
k_delay = input('k = ');
sd = input('sd = ');
end
z = menu('Change of the plant parameters', ...
'A = [1 -0.7555 0.0498]; B = [0.1 0.03]; k = 1; sd = 1; ch_time = 100; end_time = 200', ...
'Change the parameters');
disp('')
if z == 1
k_delay_n = 1;
A_n = [-0.7555 0.0498];
B_n = [zeros(1,k_delay-1), 0.1 0.03];
sd_n = 1;
ch_time = 100;
end_time = 200;
else
A_n =input('A = ');
A_n = A_n(2:length(A_n));
B_n = input('B = ');
B_n = [zeros(1,k_delay_n-1), B_n];
k_delay_n = input('k = ');
sd_n = input('sd = ');
ch_time = input('ch_time = ');
end_time = input('end_time = ');
end
nA = length(A);
nB = length(B);
nA_n = length(A_n);
nB_n = length(B_n);
% The structure and initial values of the controller parameters (polynomials R and S);
% the structure is chosen according to initial plant structure
S = [1, zeros( 1,length( A ) - 1 ) ];
R = [1, zeros( 1,length( B ) - 1 ) ];
nS = length(S);
nR = length(R);
k = k_delay;
z = menu('Controller structure and parameters', ...
'R=1; S=1; and orders according to initial plant structure', ...
'Change the parameters');
disp('')
if z == 2
S = input('S = ');
R = input('R = ');
k = input('k_delay = ');
end
randn('seed',54345);
aug1 = max(nB_n,max(nR + k,nB));
aug2 = max(nA_n,max(nS + k,nA));
u = zeros( 1, aug1 );
y = zeros( 1, aug2 );
y_ref = 0;
P = P_init * eye( nR + nS );
disp(' ');
disp('Calculations, please wait ...');
pause(0.1);
% The main loop
%______________
for t = 1:end_time
% Plant simulation
%_________________
if t < ch_time y_m = - A * y( 1:length( A ) )' + B * u( 1:length( B ) )' + sd * randn(1,1);
else y_m_n = - A_n * y( 1:length( A_n ) )' + B_n * u( 1:length( B_n ) )' + sd * randn(1,1); end;
% Recovering of vector Phi
%_________________________
Phi = [ u( k_delay:k_delay + nR-1 ), y( k_delay:k_delay + nS - 1 ) ];
% Covariance matrix updating
%___________________________
P = (1/alpha) * ( P - ( P*Phi'*Phi*P ) / ( alpha + Phi*P*Phi' ) );
% Parameters updating
%_____________________
Theta = [ R,S ] + Phi * P * ( y_m - Phi * [ R,S ]' );
R = Theta( 1:nR );
S = Theta( nR+1:nR+nS );
% Calculation of new control
%___________________________
s1 = R( 2:nR );
s2 = u( 1:nR - 1 );
if isempty( s1 ) s1 = 0; s2 = 0; end;
u_new = ( - s1 * s2' - S * [y_m, y( 1:nS - 1 ) ]' + y_ref ) / R(1);
u = [u_new , u( 1:aug1 ) ];
y = [y_m , y( 1:aug2 ) ];
% Data recording
%_______________
par_R = [ par_R,R' ];
par_S = [ par_S,S' ];
end;
% Plot of the results
%____________________
subplot(211);
plot(par_R');
title('Polynomial R parameter');
subplot(212);
plot(par_S');
title('Plynomial S parameters');
