www.gusucode.com > mpc 案例源码 matlab代码程序 > mpc/ProvidingLQRPerformanceUsingTerminalPenaltyExample.m

    %% Providing LQR Performance Using Terminal Penalty
% 
%% Define plant with one input and two outputs 

% Copyright 2015 The MathWorks, Inc.

A = [1 0;0.1 1];
B = [0.1;0.005];
C = eye(2);
D = zeros(2,1);
Ts = 0.1;
Plant = ss(A,B,C,D,Ts);
Plant.InputName = {'u'};
Plant.OutputName = {'x_1','x_2'};
%% Design LQR controller for the plant
Q = eye(2);
R = 1;
[K,Qp] = lqry(Plant,Q,R);
%% Augment the plant model
NewPlant = Plant;
cholP = chol(Qp);
set(NewPlant,'C',[C;cholP],'D',[D;zeros(2,1)],...
    'OutputName',{'x_1','x_2','Cx_1','Cx_2'})
NewPlant.InputGroup.MV = 1;
NewPlant.OutputGroup.MO = [1 2];
NewPlant.OutputGroup.UO = [3 4];
%% Create MPC controller
P = 3;
M = 3;
MPCobj = mpc(NewPlant,Ts,P,M);
%% Specify weights for MV and OV
ywt = sqrt(diag(Q))';
uwt = sqrt(diag(R))';
MPCobj.Weights.OV = [ywt 0 0];
MPCobj.Weights.MV = uwt;
MPCobj.Weights.MVrate = 1e-6;
%% Specify terminal weights
U = struct('Weight', uwt);
Y = struct('Weight', [0 0 1 1]);
setterminal(MPCobj, Y, U)
%% Remove default state estimator
setoutdist(MPCobj,'model',tf(zeros(4,1)))
setEstimator(MPCobj,[],C)
%% Compute closed-loop response with LQR 
clsys = feedback(Plant,K);

Tstop = 6;
x0 = [0.2;0.2];
[yLQR,tLQR] = initial(clsys,x0,Tstop);
%% Compute closed-loop response with MPC with terminal weights
SimOptions = mpcsimopt(MPCobj);
SimOptions.PlantInitialState = x0;
r = zeros(1,4);
[y,t,u] = sim(MPCobj,ceil(Tstop/Ts),r,SimOptions);
Cost = sum(sum(y(:,1:2)*diag(ywt).*y(:,1:2))) + sum(u*diag(uwt).*u);
%% Compute closed-loop response with standard MPC
MPCobjSTD = mpc(Plant,Ts); % Default P = 10, M = 2
MPCobjSTD.Weights.MV = uwt;
MPCobjSTD.Weights.MVrate = 1e-6;
MPCobjSTD.Weights.OV = ywt;
SimOptions = mpcsimopt(MPCobjSTD);
SimOptions.PlantInitialState = x0;
r = zeros(1,2);
[ySTD,tSTD,uSTD] = sim(MPCobjSTD,ceil(Tstop/Ts),r,SimOptions);
CostSTD = sum(sum(ySTD*diag(ywt).*ySTD)) + sum(uSTD*uwt.*uSTD);
%% Compare the responses
figure
h1 = line(tSTD,ySTD,'color','r');
h2 = line(t,y(:,1:2),'color','b');
h3 = line(tLQR,yLQR,'color','m','marker','o','linestyle','none');
xlabel('Time')
ylabel('Plant Outputs')
legend([h1(1) h2(1) h3(1)],'Standard MPC','MPC with Terminal Weights','LQR','Location','NorthEast')