%SIMEVAP Simulate evaporator under predictive control, with MPC Toolbox % % This is a script file which simulates the Newell and Lee evaporator % running under predictive control, using the MPC Toolbox. % It was used to produce the results shown in Chapter 9.2 of % `PREDICTIVE CONTROL WITH CONSTRAINTS' by J.M.Maciejowski. % % It requires the Simulink model 'evaporator.mdl' and the function % 'scmpcnl2.m', both of which are available on the book's web site. % It also requires the MPC Toolbox for Matlab. % % Edit this file to run different scenarios. % J.M. Maciejowski, 9 April 2000 % Cambridge University Engineering Dept. % Copyright (C) 2000. All Rights Reserved. %%%%% WRITE INITIAL INFO TO SCREEN: disp('Simulation of Newell and Lee Evaporator') disp(' ') disp('EVALSIM is a script file. Edit the file to change') disp(' parameters and scenarios.') disp(' ') %%%%% DEFINE VARIABLES NEEDED BY EVAPORATOR MODEL: lambda = 38.5; % kW/(kg/min) UA2 = 6.84; % kW/K Cp = 0.07; % kW/K(kg/min) lambdas = 36.6;% kW/(kg/min) rhoA = 20; % kg/m % Sample time: Ts = 1; % minute - units of time are minutes throughout. disp(['Sample time = ',num2str(Ts),' minutes']); %%%%% DEFINE INITIAL CONDITIONS FOR EVAPORATOR: % The states are ordered as follows: % (use '[sizes,x0,xstring]=evaporator' to get this information) % STATE EQUILIBRIUM VALUE % Tank Level L2 1.0 m % Composition X2 25.0 percent % Pressure P2 50.5 kPa % P100 Lag 194.7 kPa % F3 Lag 50.0 kg/min % F200 Lag 208.0 kg/min % F2 Lag 2.0 kg/min x0 = [1.0; 25.0; 50.5; 194.7; 50.0; 208.0; 2.0]; % Initial state vector %%%%% OBTAIN INITIAL LINEARISED INTERNAL MODEL: % The inputs are ordered as follows: % (ordering corresponds to Inport numbering) % INPUT EQUILIBRIUM VALUE % Control inputs: % Product flowrate F2 2.0 kg/min % Steam pressure P100 demand 194.7 kPa % Cooling water rate F200 demand 208.0 kg/min % Disturbance inputs: % Circulating rate F3 demand 50.0 kg/min % Feed flowrate F1 10.0 kg/min % Feed composition X1 5.0 percent % Feed temperature T1 40.0 deg C % Cooling water inlet temp T200 25.0 deg C % Define input values at operating point: u0 = [2.0; 194.7; 208.0; 50.0; 10.0; 5.0; 40.0; 25.0]; % Linearise: [A,B,C,D] = linmod2('evaporator',x0,u0); % Continuous-time model %Retain only 3 first outputs ( = L2, X2, P2): C = C(1:3,:); D = D(1:3,:); % Discretise model: [Ad,Bd,Cd,Dd] = ssdata(c2d(ss(A,B,C,D),Ts)); % Put model into MPC 'MOD' format: n = size(Ad,1); % Number of states nu = 3; % Number of control inputs nd = 0; % Number of measured disturances nw = 5; % Number of unmeasured disturbances nym = 3; % Number of measured outputs nyu = 0; % Number of unmeasured outputs minfo = [Ts,n,nu,nd,nw,nym,nyu]; imod = ss2mod(Ad,Bd,Cd,Dd,minfo); % Linear internal model %%%%% DEFINE CONSTRAINTS ON EVAPORATOR: % Input amplitude constraints: % (Min: 0. Max: 2*nominal operating value approx.) %ulim = [0,0,0,0,4,400,400,100]; % in MPC Toolbox format ulim = [0,0,0,4,400,400]; % Input rate constraints: Assume no constraints, since we already % have lags to represent local controllers for F2, P100, F200, F3. %ulim = [ulim,1e4,1e4,1e4,1e4]; % Set rate limits to very high values ulim = [ulim,1e4,1e4,1e4]; %Output amplitude constraints: % L2: 0 to 2m, X2: 0 to 100 percent, P2: 0 to 100 kPa. ylim = [0,0,0,2,100,100]; % in MPC Toolbox format %%%%% DEFINE PREDICTIVE CONTROL PARAMETERS: P = 30; % Prediction horizon (steps) M = 3; % Control horizon (steps) ywt = [1000,100,100]; % Weight on tracking error penalty %uwt = [1,1,1,1]; % Weight on control move penalty uwt = [1,1,1]; % Set-point trajectory: % Ramp X2 from 25 percent to 15 percent over 20 minutes, % Ramp P2 from 50.5 kPa to 70 kPa over 20 minutes, % try to hold L2 at 1.0: X2ramp = [25*ones(1,20), 25:(15-25)/20:15]; P2ramp = [50.5*ones(1,20), 50:(70-50)/20:70]; nsp = length(X2ramp); setpoints = [1.0*ones(nsp,1), X2ramp', P2ramp']; Kest = []; xm0 = []; % Default observer % No noise or measured disturbances, or unmeasured disturbances on u: z = []; d = []; wu = []; % Unmeasured disturbance: +- 10% changes on F1: F1dist = [zeros(10,1);2*ones(10,1);-ones(10,1);zeros(10,1);ones(10,1);-ones(10,1)]; nd = length(F1dist); %w = [zeros(nd,1),F1dist,zeros(nd,1),zeros(nd,1),zeros(nd,1)]; w = []; % Hard constraints: suwt = []; sywt = []; normtype = []; % Soft constraints: %suwt = 1e6; sywt = []; normtype = 'inf-norm'; disp(['Prediction horizon = ',int2str(P),' steps, ',... 'Control horizon = ',int2str(M),' steps']) disp(['Output weights : ',num2str(ywt),... '. Input weights : ',num2str(uwt)]) if isempty(normtype), disp('Hard constraints') else disp(['Soft constraints, norm = ',normtype]) end %%%%% DEFINE SIMULATION PARAMETERS: % Start and stop times: tstart = 0; tend = 100; % minutes disp(['End time = ',num2str(tend)]) %%%%% RUN SIMULATION: disp(' ') disp('To avoid re-linearisation enter interval > end time (Inf ok)') badinput = 1; while badinput, relin = input('Input interval between re-linearisations: '); if relin == 1, disp('*** Software limitation: Cannot re-linearise correctly at every step.') disp('*** Specify interval > 1') else badinput = 0; end end if relin >= tend, % No re-linearisation requested disp('Running simulation without re-linearsation') tvec = [tstart,tend]; [y,u]=scmpcnl2('evaporator',imod,ywt,uwt,M,P,tvec, ... setpoints,ulim,ylim,Kest,z,d,w,wu, ... x0,u0(1:3),[],xm0,suwt,sywt,normtype); else disp(['Re-linearising every ',num2str(relin),' minutes (simulated time)']) % Set up first segment: tvec = [tstart,relin-1]; % max() handles case relin=1 x0seg = x0; u0seg = u0(1:3); if tvec(2)+1>nsp, spseg = [setpoints(tvec(1)+1:nsp,:); setpoints(nsp,:)]; else spseg = setpoints(tvec(1)+1:tvec(2)+1,:); end y = []; u = []; while (tvec(1)+1nsp, spseg = [setpoints(tvec(1)+1:nsp,:); setpoints(nsp,:)]; else spseg = setpoints(tvec(1)+1:tvec(2)+1,:); end % Re-linearise at current state: [A,B,C,D] = linmod2('evaporator',x0seg,[u0seg;u0(4:8)]); % Continuous-time model %Retain only 3 first outputs ( = L2, X2, P2): C = C(1:3,:); D = D(1:3,:); [Ad,Bd,Cd,Dd] = ssdata(c2d(ss(A,B,C,D),Ts)); % Discrete-time imod = ss2mod(Ad,Bd,Cd,Dd,minfo); % MPC format end end %%%%% PROCESS RESULTS: tvector = (0:(size(y,1)-1))'*Ts; % Pad setpoints with final values or truncate as necessary: nt = length(tvector); if nt>nsp, lastsp = setpoints(nsp,:); setpoints = [setpoints; ones(nt-nsp,1)*lastsp]; elseif nt