Commit 4c2c8f91 authored by Patrick Göttsch's avatar Patrick Göttsch
Browse files

Gitlab Runner Update from orc-aux-files/master:Update .gitlab-ci.yml to handle...

Gitlab Runner Update from orc-aux-files/master:Update .gitlab-ci.yml to handle unzipped matlab repos for ORC
parent a1a2dd2a
function model = aircraft_model()
% ORC :: Aircraft model
%
% TUHH :: Institut for Control Systems :: Optimal and Robust Control
% Last update: 12.05.2009
% x1 = relative attitude
% x2 = forward speed
% x3 = pitch angle
% x4 = pitch rate
% x5 = vertical speed
A = [0 0 1.1320 0 -1;
0 -0.0538 -0.1712 0 0.0705;
0 0 0 1 0;
0 0.0485 0 -0.8556 -1.013;
0 -0.2909 0 1.0532 -0.6859];
% u1 = spoiler angle
% u2 = forward acceleration
% u3 = elevator angle
B = [0 0 0;
-0.12 1 0;
0 0 0;
4.419 0 -1.665;
1.575 0 -0.0732];
% y1=x1 = relative attitude
% y2=x2 = forward speed
% y3=x3 = pitch angle
C = [eye(3) zeros(3,2)];
D = zeros(3,3);
model = ss(A,B,C,D);
function GSYS = mk_LQG_gplant(model,R,Re)
% mk_LQG_gplant Make generalized plant for LQG design
%
% GSYS = mk_LQG_gplant(model,R,Re)
%
% Constructs a generalized plant GSYS for LQG design, from given model
% and weighting matrices R and Re. It is assumed that Q=C'*C, Qe=B*B'
%
% TUHH :: Institut for Control Systems :: Optimal and Robust Control
% Last update: 06.05.2009
[a,b,c,d] = ssdata(model);
n = size(a,1);
[ny,nu] = size(d);
if size(R,1) ~= nu
error('mkgplant:R','R must be an [nu x nu] matrix')
end
if size(Re,1) ~= ny
error('mkgenplant:Re','Re must be an [ny x ny] matrix')
end
A = a;
B = [zeros(n,ny) b zeros(n,ny) b]; % take Qe = b*b'
C = [c;
zeros(nu,n);
-c]; % take Q = c'*c
D = [zeros(ny+nu,ny) zeros(ny+nu,nu+ny) [zeros(ny,nu); sqrtm(R)];
eye(ny) zeros(ny,nu) sqrtm(Re) zeros(ny,nu)];
GSYS = ss(A,B,C,D);
% ORC 11.1, 14.1 :: LQG design for aircraft
%
% TUHH :: Institut for Control Systems :: Optimal and Robust Control
% Last update: 12.05.2009
%% Load the model
Gairc = aircraft_model();
[nmeas, ncont] = size(Gairc.D); % number of measured and controlled signals
%% Design parameters and generalized plant
R = eye(3); %1e-5 * mdiag(100,1,1);
% u1 = spoiler angle
% u2 = forward acceleration
% u3 = elevator angle
Re = eye(3);
% take Q = c'*c and Qe = b*b'
% Try creating your own function mk_LQG_gplant()
GSYS = mk_LQG_gplant( Gairc, R, Re );
%% Controller design
K = h2syn( GSYS, nmeas, ncont )
% Note: The old Matlab "Robust control toolbox" requires ltisys type of
% systems. If you haven't configured your Matlab to use the new toolbox you
% would need the function mk_LQG_old to design your controller
% K = mk_LQG_old(GSYS,[nmeas+ncont nmeas ncont]);
%% Simulation
% Using Matlab - fast
Gyr = feedback ( Gairc*K, eye(nmeas) );
Gur = feedback ( K, Gairc );
[y,t] = step(Gyr, 10);
[u,t] = step(Gur, t);
% Select the first input channel: step to spoiler angle change
y = y(:,:,1);
u = u(:,:,1);
% % Or with Simulink
% rswitch = 1; % reference
% dswitch = 0; % disturbance to be used later
% nswitch = 0; % noise
%
% sim('aircraft_sim',10); % step to spoiler angle change
% t = output.time;
% y = output.signals.values;
% u = input.signals.values;
% Plot
figure
subplot(211); plot(t,y); title('System outputs'); grid
subplot(212); plot(t,u); title('Control inputs') , legend({'u1','u2','u3'}); grid
%% Sigma plot for ORC 14.1
figure
sigma(Gyr)
% ORC 15.2 :: LQG design for aircraft
%
% TUHH :: Institut for Control Systems :: Optimal and Robust Control
% Last update: 20.05.2009
% Note: This is as in Exercise 11.1, but without the reference signal
%% Load the model
Gairc = aircraft_model();
[A,B,C,D] = ssdata(Gairc);
[nmeas, ncont] = size(D); % number of measured and controlled signals
nstates = size(A,1);
%% Generalized plant
% [ A | Bw B ]
% GSYS = [ ----+----------]
% [ Cz | 0 Dzu ]
% [ -C | Dvw 0 ]
% Design 1
Cz = 2*[C; zeros(3,nstates)];
Dzu = [zeros(3); eye(3)];
Bw = 2*[B zeros(nstates,3)];
Dvw = [zeros(3) eye(3)];
GSYS1 = ss( A, [Bw B], [Cz; -C], [zeros(6) Dzu; Dvw zeros(3)]);
% Design 2
Cz = 50*[C; zeros(3,nstates)];
Dzu = [zeros(3); eye(3)];
Bw = 50*[B zeros(nstates,3)];
Dvw = [zeros(3) eye(3)];
GSYS2 = ss( A, [Bw B], [Cz; -C], [zeros(6) Dzu; Dvw zeros(3)]);
%% Controller design
[K1,CL1,h1] = h2syn( GSYS1, nmeas, ncont );
[K2,CL2,h2] = h2syn( GSYS2, nmeas, ncont );
fprintf('H2 norm with K1 = %f\nH2 norm with K2 = %f\n',h1,h2)
%% Singular value plots
Gyr1 = feedback(Gairc*K1,eye(3));
Gyr2 = feedback(Gairc*K2,eye(3));
save orc15_design2 K2 Gairc CL2 Gyr2 % save for later
figure(1)
sigma(Gyr1,'r--',Gyr2,'b-')
legend('Design 1','Design 2')
%% Input directions
SteadyStatefrsp = freqresp(Gyr1, 0);
[Uss, Sss, Vss] = svd(SteadyStatefrsp)
% Vss = Vss;
Tsim = 12;
t = 0:0.01:Tsim;
% Select the columns of Vss that correspond to the singular value closest
% to 1 (best input direction) or farthest from 1 (worst input direction)
rss_worst = Vss(:,1)*ones(1, length(t));
rss_best = Vss(:,2)*ones(1, length(t));
[yss_worst, t] = lsim(Gyr1, rss_worst, t);
[yss_best , t] = lsim(Gyr1, rss_best , t);
% Plot
figure(2)
subplot(211); plot(t, yss_worst,'r', t, rss_worst, 'k'); title('System outputs, worst input direction'); grid
subplot(212); plot(t, yss_best ,'r', t, rss_best , 'k'); title('System outputs, best input direction'); grid
%% Simulation
Tsim = 8;
% With Matlab
Gur1 = feedback ( K1, Gairc );
Gur2 = feedback ( K2, Gairc );
[y1,t1] = step(Gyr1, Tsim);
[u1,t1] = step(Gur1, t1);
[y2,t2] = step(Gyr2, Tsim);
[u2,t2] = step(Gur2, t2);
% Select the first input channel: step to spoiler angle change
y1 = y1(:,:,1); u1 = u1(:,:,1);
y2 = y2(:,:,1); u2 = u2(:,:,1);
% Plot
figure(3)
subplot(211); plot(t1,y1,'r--',t2,y2,'b-'); title('System outputs'); grid
subplot(212); plot(t1,u1,'r--',t2,u2,'b-'); title('Control inputs'); grid
xlabel('Design 1 [red --], Design 2 [blue -]')
%% Or with Simulink
rswitch = 0; % reference
dswitch = 1; % disturbance to be used later
nswitch = 1; % noise
K = K1;
sim('aircraft_sim', Tsim); % step to spoiler angle change
t1 = output.time;
y1 = output.signals.values;
u1 = input.signals.values;
K = K2;
sim('aircraft_sim', Tsim); % step to spoiler angle change
t2 = output.time;
y2 = output.signals.values;
u2 = input.signals.values;
% Plot
figure(3)
subplot(211); plot(t1,y1,'r--',t2,y2,'b-'); title('System outputs'); grid
subplot(212); plot(t1,u1,'r--',t2,u2,'b-'); title('Control inputs'); grid
xlabel('Design 1 [red --], Design 2 [blue -]')
% ORC 16.2 :: Sensitivity
%
% TUHH :: Institut for Control Systems :: Optimal and Robust Control
% Last update: 28.05.2009
%% Repeat the LQG controller design for the aircraft:
orc11_1_and_14_1_LQG
close all
% now we have: K, Gairc
%% Plot the singular value plot for the sensitivity
I = eye(3);
T = feedback ( Gairc*K, I );
S = feedback( I, Gairc*K );
figure(1), sigma(T,'b-')
figure(2), sigma(S, 'r-')
%% Find the best and worst input directions for low frequencies
w0 = 0.001;
S0 = freqresp(S,w0); % the same as C*inv(j*w0-A)*B+D
[u,s,v] = svd(S0)
%%
disp('Best direction')
v(:,3)
disp('System response to best direction')
S0 * v(:,3)
u(:,3)*s(3,3)
disp(' ')
disp('Worst direction')
v(:,1)
disp('System response to worst direction')
S0 * v(:,1)
%% Several simulations to prove the above results
% Note! the best direction is [1; 0; 0], i.e. first input
figure(3)
step( S(:,1),'r-', S(:,2),'g-', S(:,3),'g-',25);
legend('Best direction','Other directions')
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment