%% State Space MPC Tutorial
% This document explains how to use the setup function and online
% controller returned by ssmpcsetup.
%% State-space MPC set-up
% The online controller has to be set-up before use. To set-up the state
% space MPC controller, the user has to supply the state-space model,
% represented by A, B, C and D matrices (Note: it should be in
% discrete-time), the predication horizon, P, moving horizon, M and
% performance weights, Q and R. The default initial state and input are set
% to zero.
%
% SSMPC=MPCSETUP(A,B,C,D,P,M,Q,R,X0,U0);
%% Online controller
% The returned function handle from the MPC setup program is an online
% controller, SSMPC. The controller is called by proving two inputs:
% current measurement, Y and future reference, Ref. On return, it produces
% the optimal input, U for next step:
%
% U = SSMPC(Y,Ref);
sys='OR33';
Ts=30;
[A,B,C,D]=dlinmod(sys,Ts);
%% MPC parameters
% The MPC controller is configured with following parameters.
% Prediction horizon and moving horizon
p=112;
m=70;
% Performance wights
Q=1.5*eye(3*p);
R=eye(3*m);
%% MPC set-up
% The MPC controller is set-up by calling SSMPCSETUP:
ssmpc=mpcsetup(A,B,C,D,p,m,Q,R);
%% Simulation
% 150 seconds (1500 sampling intervals) simulation is conducted with
% several setpoint changes and random cooling water temperature changes
% within positive and negative 1 degree.
% Simulation length and variables for results
N=250;
x0=zeros(10,1);
Y=zeros(N,3);
U=zeros(N,3);
% Predefined reference
T=zeros(N,3);
t1=50;
t2=100;
t3=150;
t4=200;
t5=N;
T(t1:N,:)=1;
T(t2:N,:)=2;
T(t3:N,:)=3;
T(t4:N,:)=4;
% Simulation
for k=1:N
% Process disturbances
w=0; %Bd*(rand(1,1)-0.5)*2;
% Measurements noise
v=0.01*randn(3,1);
% actual measurement
y=C*x0+v;
% online controller
u=ssmpc(y,T(k:end,:));
% plant update
x0=A*x0+B*u+w;
% save results
Y(k,:)=y';
U(k,:)=u';
end
%% Results
% The simulation results are summarized in two sub-plots.
t=(0:N-1);
figure(1)
subplot(311)
plot(t,Y(:,1),t,T(:,1),'r--','linewidth',2)
title('output and setpoint')
ylabel('ETHANOL mol.fract in the D, %')
legend('y_1','Ref','location','southeast')
subplot(312)
plot(t,Y(:,2),t,T(:,2),'r--','linewidth',2)
title('output and setpoint')
ylabel('ETHANOL mol.fract in the SS, %')
legend('y_2','Ref','location','southeast')
subplot(313)
plot(t,Y(:,3),t,T(:,3),'r--','linewidth',2)
title('output and setpoint')
ylabel('TEMP in the TRAY 19, C')
legend('y_3','Ref','location','southeast')
figure(2)
subplot(311)
stairs(t,U(:,1),'linewidth',2)
legend('u_1','location','southeast')
title('input')
ylabel('reflux flowrate, m^3/s')
xlabel('time, s')
subplot(312)
stairs(t,U(:,2),'linewidth',2)
legend('u_2','location','southeast')
title('input')
ylabel('sidestream product flowrate')
xlabel('time, s')
subplot(313)
stairs(t,U(:,3),'linewidth',2)
legend('u_3','location','southeast')
title('input')
ylabel('reboiler stream pressure')
xlabel('time, s')
n=length(t);
% ISE (integral square error calculation for y1
ISE=0;
IAE=0;
for i=2:t1
ISE(i)=ISE(i-1)+(0-Y(i-1,1))^2;
IAE(i)=IAE(i-1)+abs(0-Y(i-1,1));
end
ISE0y1=ISE(i);
IAE0y1=IAE(i);
ISE(t1)=0;
IAE(t1)=0;
for i=t1+1:t2
ISE(i)=ISE(i-1)+(T(i-1)-Y(i-1,1))^2;
IAE(i)=IAE(i-1)+abs(T(i-1)-Y(i-1,1));
end
ISEt1y1=ISE(i);
ISE(t1)=0;
IAEt1y1=IAE(i);
IAE(t1)=0;
for i=t2+1:t3
ISE(i)=ISE(i-1)+(T(i-1)-Y(i-1,1))^2;
IAE(i)=IAE(i-1)+abs(T(i-1)-Y(i-1,1));
end
ISEt2y1=ISE(i);
ISE(t2)=0;
IAEt2y1=IAE(i);
IAE(t2)=0;
for i=t3+1:t4
ISE(i)=ISE(i-1)+(T(i-1)-Y(i-1,1))^2;
IAE(i)=IAE(i-1)+abs(T(i-1)-Y(i-1,1));
end
ISEt3y1=ISE(i);
ISE(t3)=0;
IAEt3y1=IAE(i);
IAE(t3)=0;
for i=t4+1:t5
ISE(i)=ISE(i-1)+(T(i-1)-Y(i-1,1))^2;
IAE(i)=IAE(i-1)+abs(T(i-1)-Y(i-1,1));
end
ISEt4y1=ISE(i);
ISEfory1=[ISE0y1 ISEt1y1 ISEt2y1 ISEt3y1 ISEt4y1]
IAEt4y1=IAE(i);
IAEfory1=[IAE0y1 IAEt1y1 IAEt2y1 IAEt3y1 IAEt4y1]
% ISE (integral square error calculation for y2
ISE=0;
IAE=0;
for i=2:t1
ISE(i)=ISE(i-1)+(0-Y(i-1,2))^2;
IAE(i)=IAE(i-1)+abs(0-Y(i-1,2));
end
ISE0y2=ISE(i);
IAE0y2=IAE(i);
ISE(t1)=0;
IAE(t1)=0;
for i=t1+1:t2
ISE(i)=ISE(i-1)+(T(i-1)-Y(i-1,2))^2;
IAE(i)=IAE(i-1)+abs(T(i-1)-Y(i-1,2));
end
ISEt1y2=ISE(i);
ISE(t1)=0;
IAEt1y2=IAE(i);
IAE(t1)=0;
for i=t2+1:t3
ISE(i)=ISE(i-1)+(T(i-1)-Y(i-1,2))^2;
IAE(i)=IAE(i-1)+abs(T(i-1)-Y(i-1,2));
end
ISEt2y2=ISE(i);
ISE(t2)=0;
IAEt2y2=IAE(i);
IAE(t2)=0;
for i=t3+1:t4
ISE(i)=ISE(i-1)+(T(i-1)-Y(i-1,2))^2;
IAE(i)=IAE(i-1)+abs(T(i-1)-Y(i-1,2));
end
ISEt3y2=ISE(i);
ISE(t3)=0;
IAEt3y2=IAE(i);
IAE(t3)=0;
for i=t4+1:t5
ISE(i)=ISE(i-1)+(T(i-1)-Y(i-1,2))^2;
IAE(i)=IAE(i-1)+abs(T(i-1)-Y(i-1,2));
end
ISEt4y2=ISE(i);
ISEfory2=[ISE0y2 ISEt1y2 ISEt2y2 ISEt3y2 ISEt4y2]
IAEt4y2=IAE(i);
IAEfory2=[IAE0y2 IAEt1y2 IAEt2y2 IAEt3y2 IAEt4y2]
% ISE (integral square error calculation for y3
ISE=0;
IAE=0;
for i=2:t1
ISE(i)=ISE(i-1)+(0-Y(i-1,3))^2;
IAE(i)=IAE(i-1)+abs(0-Y(i-1,3));
end
ISE0y3=ISE(i);
IAE0y3=IAE(i);
ISE(t1)=0;
IAE(t1)=0;
for i=t1+1:t2
ISE(i)=ISE(i-1)+(T(i-1)-Y(i-1,3))^2;
IAE(i)=IAE(i-1)+abs(T(i-1)-Y(i-1,3));
end
ISEt1y3=ISE(i);
ISE(t1)=0;
IAEt1y3=IAE(i);
IAE(t1)=0;
for i=t2+1:t3
ISE(i)=ISE(i-1)+(T(i-1)-Y(i-1,3))^2;
IAE(i)=IAE(i-1)+abs(T(i-1)-Y(i-1,3));
end
ISEt2y3=ISE(i);
ISE(t2)=0;
IAEt2y3=IAE(i);
IAE(t2)=0;
for i=t3+1:t4
ISE(i)=ISE(i-1)+(T(i-1)-Y(i-1,3))^2;
IAE(i)=IAE(i-1)+abs(T(i-1)-Y(i-1,3));
end
ISEt3y3=ISE(i);
ISE(t3)=0;
IAEt3y3=IAE(i);
IAE(t3)=0;
for i=t4+1:t5
ISE(i)=ISE(i-1)+(T(i-1)-Y(i-1,3))^2;
IAE(i)=IAE(i-1)+abs(T(i-1)-Y(i-1,3));
end
ISEt4y3=ISE(i);
ISEfory3=[ISE0y3 ISEt1y3 ISEt2y3 ISEt3y3 ISEt4y3]
IAEt4y3=IAE(i);
IAEfory3=[IAE0y3 IAEt1y3 IAEt2y3 IAEt3y3 IAEt4y3]
ISEfory1TOTAL=ISE0y1+ISEt1y1+ISEt2y1+ISEt3y1+ISEt4y1
ISEfory2TOTAL=ISE0y2+ISEt1y2+ISEt2y2+ISEt3y2+ISEt4y2
ISEfory3TOTAL=ISE0y3+ISEt1y3+ISEt2y3+ISEt3y3+ISEt4y3
mpc.rar_DMC_MPC_MPC fuzzy_mpc tuning_pon
版权申诉
187 浏览量
2022-07-14
13:38:44
上传
评论
收藏 9KB RAR 举报
JonSco
- 粉丝: 69
- 资源: 1万+