Commit 060d6790 authored by Prima Aditya's avatar Prima Aditya
Browse files

Gitlab Runner Update from cstd-aux-files/master:Merge branch...

Gitlab Runner Update from cstd-aux-files/master:Merge branch 'Changes_according_to_new_script' into 'master'
parent 9c9500a1
%% Controller design for process evaporator
% This skript uses a model of the system, obtained with the identification
% toolbox. The model should be named mod1
Ts = 0.05; % sampling time
[Ad,Bd,Cd,Dd] = ssdata(mod1);
SysD = ss(Ad,Bd,Cd,Dd, Ts); % discrete model
% This re-definition is needed, because of the special idss-form of mod1
n = size(Ad,1) % order of the model
[no, ni] = size(Dd) % number of inputs and outputs
%% Define the continuous time model
SysCns=d2c(SysD,'zoh');
%introduce scaling here if required
umax=[1 1 1]; % See Control objective 2.c
ymax=[.1 2 2]; % See Control objective 2.a [0.1 2 2];
scale_inc =diag(umax);
scale_outc=diag(1./ymax);
SysC = scale_outc * SysCns * scale_inc; % scaled continuous time model
[A,B,C,D] = ssdata(SysC);
%% Design a state-feedback and obsever
% FIXME FIXME FIXME FIXME
% This is already implemented in
controllerdesign2
% but as a self-test try writing by yourself the files for controller and
% observer design
% The script should return
% f0, fI, Aobs, Bobs
%% Obtain model of the closed-loop system with the linear model
Tmax = 75; % max simulation time (in minutes)
[Acllqg,Bcllqg,Ccllqg,Dcllqg] = linmod('evaplinmod');
Scllqg = ss(Acllqg,Bcllqg,Ccllqg,Dcllqg);
%% Simultate the closed-loop system
T1 = (0:Ts:Tmax)';
[Yr1q,T1]=step(Scllqg(:,1),T1);
[Yr2q,T1]=step(Scllqg(:,2),T1);
[Yr3q,T1]=step(Scllqg(:,3),T1);
%% Calculate initial states for the non-linear model
% The non-linear model should always start from a stationary point
% calculate steady state for observer
xobs0 = -inv(Aobs)*Bobs*[inv(scale_inc)*u0; scale_outc*y0];
% xobs0 = -inv(Aobs)*Bobs*[u0; y0];
% steady state for integrator
xint0 = inv(fI)*(inv(scale_inc)*u0 - f0*xobs0);
%% Simulate the non-linear model and plot the outputs and control inputs of
%% both the non-linear and the linear model
% step_time = 1; % step after ... minutes
%
% comp_input = 'level'; % Compensate level
% simandplot
%
%
% comp_input = 'pressure';
% simandplot
% %
% comp_input = 'temp';
% simandplot
%
step_time = 20;
comp_input = 'dist';
simandplot
% get range of movement for inputs and outputs % outputs: 1:3, inputs 4:6
rangevec = range(outputdat)
stdvec = std(outputdat(end-200:end,:)) % standard deviation outputs/inputs
%% simple function, finding the index of the settling time moment
settime = @(vec,dev) find(abs(vec-vec(1))>(0.05+2*dev), 1, 'last');
stvec_level = T1( settime(outputdat(:,1),1*stdvec(1)) ) - step_time
stvec_press = T1( settime(outputdat(:,2),1*stdvec(2)) ) - step_time
stvec_temp = T1( settime(outputdat(:,3),1*stdvec(3)) ) - step_time
This diff is collapsed.
%% LQG controller design
% Define augmented plant
Aa = [A zeros(n,ni);
-C zeros(ni)];
Ba = [B;
-D];
% Define the weighting matrices - for states and integrators
% Alternative to defining R, Q, RN and QN, just load a setting from
tuningparset
% EAch setting is defined with different color, so one can compare them
% Design state-feedback
F = lqr(Aa,Ba,Q,R);
f0 = -F(:,1:n); % states-gain
fI = -F(:,n+(1:ni)); % integrator-gain
Lt = lqr(A',C',QN,RN);
L = -Lt';
% Build the observer matrices
Aobs = A + L*C;
Bobs = [B -L];
% discretise the dynamic controller
% get controller: 6 inputs (3 setpoints, 3 measurements)
% 3 outputs: inputs to plant
Ts = 0.05
[Ac,Bc,Cc,Dc]=linmod('controller');
% initialise values
% u1-u3: setpoints=y0
% u4-u6: measured values=y0
% y=u0
u0cont = [scale_outc*y0;scale_outc*y0];
y0cont = inv(scale_inc)*u0;
xcont0 = pinv([Ac;Cc])*[-Bc*u0cont; (y0cont-Dc*u0cont)];
Ccont = ss(Ac,Bc,Cc,Dc);
[Cdisc,M] = c2d(Ccont,Ts,'zoh');
% calculate initial state of discrete time state space model
dxcont0 = M*[xcont0; u0cont];
%% Run model to get graphs of response to disturbance
step_time = 20;
ns = length(T1); % number of samples to simulate
sp1 = ones(ns,1)*scale_outc(1,1)*y0(1);
sp2 = ones(ns,1)*scale_outc(2,2)*y0(2);
sp3 = ones(ns,1)*scale_outc(3,3)*y0(3);
dist = 2*ones(ns,1);
dist(T1>step_time) = 2.4;
% Set up the input vector data structure to the model
Usig.time = T1;
% Usig.signals.values = [sp1 sp2 sp3 dist];
% Usig.signals.dimensions = 4;
Usig.signals(1).values = SP(:,1);
Usig.signals(1).dimensions = 1;
Usig.signals(2).values = SP(:,2);
Usig.signals(2).dimensions = 1;
Usig.signals(3).values = SP(:,3);
Usig.signals(3).dimensions = 1;
Usig.signals(4).values = dist;
Usig.signals(4).dimensions = 1;
% Simulate the closed-loop with the discrete model
sim('evapmodd',T1,[],Usig);
figure(9); % outputs - scale back the non-linear ones
subplot(311), plot(T1,outputdat(:,1) - scale_outc(1,1)*y0(1)), title(['Comparison ' comp_input ' outputs'])
hold on, ylabel('level')
subplot(312), plot(T1,outputdat(:,2) - scale_outc(2,2)*y0(2))
hold on, ylabel('pressure')
subplot(313), plot(T1,outputdat(:,3) - scale_outc(3,3)*y0(3))
hold on, ylabel('temp')
figure(10) % inputs - scale back the non-linear ones
subplot(311), plot(T1,outputdat(:,4) - 1/scale_inc(1,1)*u0(1)), title(['Comparison ' comp_input ' inputs'])
hold on, ylabel('fluid out')
subplot(312), plot(T1,outputdat(:,5) - 1/scale_inc(2,2)*u0(2))
hold on, ylabel('gas out')
subplot(313), plot(T1,outputdat(:,6) - 1/scale_inc(3,3)*u0(3))
hold on, ylabel('heat fuild')
function [sys,x0,str,ts] =evap2s(t,x,u,flag,x0evap,pvec)
% This is a S-function, modelling the evaporator
% It implements the equations of evaporator and is used by Simulink during
% the simulation of the non-linear model
switch flag,
%%%%%%%%%%%%%%%%%%
% Initialization %
%%%%%%%%%%%%%%%%%%
case 0,
[sys,x0,str,ts]=mdlInitializeSizes(x0evap);
%%%%%%%%%%%%%%%
% Derivatives %
%%%%%%%%%%%%%%%
case 1,
sys=mdlDerivatives(t,x,u,pvec);
%%%%%%%%%%%
% Outputs %
%%%%%%%%%%%
case 3,
sys=mdlOutputs(t,x,u,pvec);
%%%%%%%%%%%%%%%%%%%
% Unhandled flags %
%%%%%%%%%%%%%%%%%%%
case { 2, 4, 9 },
sys = [];
%%%%%%%%%%%%%%%%%%%%
% Unexpected flags %
%%%%%%%%%%%%%%%%%%%%
otherwise
error(['Unhandled flag = ',num2str(flag)]);
end
% end csfunc
%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts]=mdlInitializeSizes(x0evap);
sizes = simsizes;
sizes.NumContStates = 4;
sizes.NumDiscStates = 0;
sizes.NumOutputs = 3;
sizes.NumInputs = 5;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1;
sys = simsizes(sizes);
x0 =x0evap;
str = [];
ts = [0 0];
% end mdlInitializeSizes
%
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%
function sys=mdlDerivatives(t,x,u,pvec)
sys=evapderiv(t,x,u,pvec);
% end mdlDerivatives
%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
%
function sys=mdlOutputs(t,x,u,pvec)
cj=pvec(7);
vj=pvec(3);
%output scaling
dLmax=0.25;
dPmax=0.1;
dTmax=0.01;
%output 1 : normalised volume (scaled)
sys(1)=x(1)/dLmax;
%output 2:normalised pressure (scaled)
sys(2)=x(2)*x(3)/x(1)*1/(2-x(1))/dPmax;
%output 3: normalised temperature (scaled)
sys(3)=x(3)/x(1)/dTmax;
% end mdlOutputs
function deriv=evapderiv(t,x,u,pvec)
% Computes the derivatives needed for the S-function, implementing the
% nonlinear equations of the evaporator
% pvec = [taul,taug,tauj,UA,p1,kvap,cj,Tjin,lat,hin,t1,dH];
taul=pvec(1);
taug=pvec(2);
vj=pvec(3);
UA=pvec(4);
pnorm=pvec(5);
kvap=pvec(6);
cj=pvec(7);
Tjin=pvec(8);
lat=pvec(9);
%hin=pvec(10);
tnorm=pvec(11);
dH=pvec(12);
%input scaling
dflmax=0.5;
dfgmax=0.5;
dfjmax=0.03;
scalin=diag([dflmax,dfgmax,dfjmax]);
fvec=scalin*u(2:4);
fin=u(1)*dflmax;
flout=fvec(1);
fgout=fvec(2);
fjin=fvec(3);
hin=u(5);
%x2 is PV/Tl
vl=x(1);
H=x(3);
pvgotl=x(2);
Tl=H/vl;
%liquid volume and gas volume are identical
vg=2-vl;
p=x(2)*x(3)/x(1)*1/vg;
Pvap=pnorm*exp(dH*(1/tnorm-1/Tl));
flg=kvap*fin*(Pvap-p);
hlout=Tl;
Tj=x(4)/(cj*vj);
hjout=Tj*cj;
hjin=Tjin*cj;
Qjv=UA*(Tj-Tl);
deriv(1)=1/taul*(fin-flout-flg);
deriv(2)=1/taug*(flg-fgout);
deriv(3)=1/taul*(fin*hin-flout*hlout-lat*flg+Qjv);
deriv(4)=1/taul*(fjin*(hjin-hjout)-Qjv);
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
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