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')