www.pudn.com > nnctrl_v5.zip > lintest.m, change:1997-06-06,size:8735b

```% PROGRAM FOR DEMONTRATING CONTROL BASED ON INSTANTANEOUS LINEARIZATION
%
% Programmed by Magnus Norgaard, IAU/IMM, Technical Univ. of Denmark
% LastEditDate: Feb. 20, 1996
close all
StopDemo=0;
figure
guihand=gcf;
for k=1:1, %dummy loop

% >>>>>>>>>>>>>>>>  BUILD GUI INTERFACE  <<<<<<<<<<<<<<<<<
[guihand,edmulti,contbut,quitbut]=pmnshow;
set(guihand,'Name','Control based on instantaneous linearization');

% >>>>>>>>>>>>>>>>  SCREEN 1  <<<<<<<<<<<<<<<<<
s0='1';
s1='The purpose of this demo is to show how the instantaneous';
s2='linearization principle can be used in the design of';
s3='control systems for nonlinear processes. Two types of';
s4='designs will be investigated: pole placement with zero';
s5='cancellation and pole placement without zero cancellation.';
s6='The process in question is a spring-mass-damper system';
s7='with a hardening spring: y"(t) + y''(t) + y(t) + y(t)^{3} = u(t)';
smat=str2mat(s0,s1,s2,s3,s4,s5,s6,s7);
pmnshow(smat,guihand,edmulti,contbut,quitbut);
if StopDemo==1, close all, break; end

% >>>>>>>>>>>>>>>>  SCREEN 2  <<<<<<<<<<<<<<<<<
% -- Generate data --
N2=length(U);
N1=floor(N2/2);
Y1 = Y(1:N1)';
U1 = U(1:N1)';
Y2 = Y(N1+1:N2)';
U2 = U(N1+1:N2)';
s0='2';
s1='Before we can apply the controller design we need a neural';
s2='network model of the process. To create this we must make';
s3='an experiment and collect a set of data describing the';
s4='process over its entire range of operation. Such an';
s5='experiment has been simulated in advance with the function';
s6='"experim." The plot above shows the data set.';
smat=str2mat(s0,s1,s2,s3,s4,s5,s6);

subplot(411)
plot(U1); grid
axis([0 N1 min(U1) max(U1)])
title('Input and output sequence')
subplot(412)
plot(Y1); grid
axis([0 N1 min(Y1) max(Y1)])
xlabel('time (samples)')
drawnow
pmnshow(smat,guihand,edmulti,contbut,quitbut);
if StopDemo==1, close all, break; end

% >>>>>>>>>>>>>>>>  SCREEN 3  <<<<<<<<<<<<<<<<<
s0='3';
s1='To identify the neural network model we will use the';
s2='function "nnarx" from the NNSYSID-toolbox. Since it''s';
s3='a second order process we will use as regressors two';
s4='past outputs and two past controls. Furthermore we choose';
s5='a network architecture with five hidden "tanh" units and';
s6='one linear output.';
subplot(411);delete(gca);subplot(412);delete(gca)
subplot('position',[0.1 0.55 0.45 0.38]);
drawnet(ones(7,5),ones(1,8),eps,['y(t-1)';'y(t-2)';'u(t-1)';'u(t-2)'],'yhat(t)');
title('Network architecture')
smat=str2mat(s0,s1,s2,s3,s4,s5,s6);
pmnshow(smat,guihand,edmulti,contbut,quitbut);
if StopDemo==1, close all, break; end

% >>>>>>>>>>>>>>>>  SCREEN 4  <<<<<<<<<<<<<<<<<
% ----- Train network -----
s0='4';
s1=[];
s2='    >> Training process in action!! <<';
s3=[];
s4=[];
s5='We run up to 200 iterations so you may have to';
s6='wait for a while.';
smat=str2mat(s0,s1,s2,s3,s4,s5,s6);
set(edmulti,'String',smat);
drawnow
trparms = [200 0 1 0];
NetDef  = ['HHHHH';'L----'];
NN=[2 2 1];
[W1,W2]=nnarx(NetDef,NN,[],[],trparms,Y1,U1);
save forward2 W1 W2 NetDef NN
delete(gca);
subplot('position',[0.1 0.55 0.45 0.38]);
drawnet(W1,W2,eps,['y(t-1)';'y(t-2)';'u(t-1)';'u(t-2)'],'yhat(t)');
title('Trained network')
if StopDemo==1, close all, break; end

% >>>>>>>>>>>>>>>>  SCREEN 5  <<<<<<<<<<<<<<<<<
s0='5';
s1='The network has now been trained and we are ready to';
s2='simulate the control system. First we will try a';
s3='pole placement strategy where we cancel the zero.';
s4='Let''s choose the desired closed loop system as:';
s5='                0.09z^{3}';
s6='H(z)=----------------------------';
s7='        z^{3} - 1.4z^{2} + 0.49z';
s8='corresponding to two poles in z=0.7 and one in z=0.';
smat=str2mat(s0,s1,s2,s3,s4,s5,s6,s7,s8);
pmnshow(smat,guihand,edmulti,contbut,quitbut);
if StopDemo==1, close all, break; end

% >>>>>>>>>>>>>>>>  SCREEN 6  <<<<<<<<<<<<<<<<<
figure('Units','Centimeters','Position',[1.5 1.5 10 1.5]);
pp=1;
lincon
close
subplot(411)
plot([0:samples-1],[ref_data y_data ym_data]); grid
axis([0 samples -2 2])
title('Reference, output and desired output')
subplot(412)
plot([0:samples-1],u_data);
axis([0 samples min(u_data) max(u_data)]); grid
title('Control signal')
xlabel('time (samples)')
drawnow
s0='6';
s1='Very nice model-following is clearly achieved. However';
s2='the control signal is very active because the zero we';
s3='cancel is close to -1. To avoid this we instead';
s4='choose not to let the controller cancel the zero.';
s5='Let''s use the same desired model as before but introduce';
s6='a dead-beat observer polynomial for causality reasons:';
s7='Ao(z)=z';
smat=str2mat(s0,s1,s2,s3,s4,s5,s6,s7);
pmnshow(smat,guihand,edmulti,contbut,quitbut);
if StopDemo==1, close all, break; end

% >>>>>>>>>>>>>>>>  SCREEN 7  <<<<<<<<<<<<<<<<<
figure('Units','Centimeters','Position',[1.5 1.5 10 1.5]);
pp=2;
lincon
close
subplot(411)
plot([0:samples-1],[ref_data y_data]); grid
axis([0 samples -2 2])
title('Reference and output signal')
subplot(412)
plot([0:samples-1],u_data);
axis([0 samples min(u_data) max(u_data)]); grid
title('Control signal')
xlabel('time (samples)')
drawnow
s0='7';
s1='It is seen that the response looks more or less as';
s2='before but that the control signal now is much more';
s3='smooth.';
smat=str2mat(s0,s1,s2,s3);
pmnshow(smat,guihand,edmulti,contbut,quitbut);
if StopDemo==1, close all, break; end

% >>>>>>>>>>>>>>>>  SCREEN 8  <<<<<<<<<<<<<<<<<
subplot(411)
plot([0:samples-1],B_data); grid
axis([0 samples min(min(B_data(3:samples,:))) max(max(B_data))])
title('Numerator (above) and denominator (below) coefficients')
subplot(412)
plot([0:samples-1],A_data(:,2:3));
axis([0 samples min(min(A_data(3:samples,2:3))) max(max(A_data))]); grid
title('Control signal')
xlabel('time (samples)')
drawnow
s0='8';
s1='A very nice feature of instantaneous linearization is';
s2='that it provides an excellent understanding of the';
s3='dynamics of the process. As shown above we can, for';
s4='example, plot the coefficients of the extracted linear';
s5='models.';
s6='   Another possibility is to plotpoles and zeros in';
s7='the complex plane.';
smat=str2mat(s0,s1,s2,s3,s4,s5,s6,s7);
pmnshow(smat,guihand,edmulti,contbut,quitbut);
if StopDemo==1, close all, break; end

% >>>>>>>>>>>>>>>>  SCREEN 9  <<<<<<<<<<<<<<<<<
subplot(411);delete(gca);subplot(412);delete(gca)
subplot('position',[0.1 0.55 0.45 0.38]);
Zmat = zeros(samples-2,nb-1);
for k=3:samples,
Zmat(k-2,:)=roots(B_data(k,:))';
end
ZPmat = zeros(samples-2,na);
for k=3:samples,
ZPmat(k-2,:)=roots(A_data(k,:))';
end
plot(ZPmat,'x');
title('Poles (x) and zeros (o) of extracted linear models');grid
axis([-1 1 -1 1])
axis('equal')
hold
plot(Zmat,Zmat*0,'o');
t=-pi:0.05:pi;
plot(sin(t),cos(t))
hold off
drawnow
s0='9';
s1='From this plot it is quite clear why there were problems';
s2='when the controller canceled the zero. Obviously it is';
s3='very close to the unity circle. We can also see that the';
s4='poles varies quite a lot. To get a better understanding';
s5='of how the pole locations affect the behavior of the';
s6='process we can calculate damping factor and natural';
s7='frequency.';
smat=str2mat(s0,s1,s2,s3,s4,s5,s6,s7);
pmnshow(smat,guihand,edmulti,contbut,quitbut);
if StopDemo==1, close all, break; end

% >>>>>>>>>>>>>>>>  SCREEN 10  <<<<<<<<<<<<<<<<<
Cmat = zeros(samples-2,na-1);
for k=3:samples,
r=roots(A_data(k,:))';
[n,m]=sort(real(r));
[MAG,a1,a2]=ddamp(r(m)',0.2);
Cmat(k-2,1)=a1(1);
Cmat(k-2,2)=a2(1);
end
Cmat=real(Cmat);
subplot(411)
plot(2:samples-1,Cmat(:,1));grid
title('Natural frequency (rad/s) and damping factor')
axis([0 samples min(Cmat(:,1)) max(Cmat(:,1))])
subplot(412)
plot(2:samples-1,Cmat(:,2));grid
axis([0 samples 0 1])
xlabel('time (samples)')
drawnow
s0='10';
s1='Obviously the process becomes less damped and the';
s2='natural frequency larger when the magnitude of';
s3='the output increases.';
s4=[];
s5=[];
s6='              >>  THE END <<';
smat=str2mat(s0,s1,s2,s3,s4,s5,s6);
set(edmulti,'String',smat);
drawnow
end

```