mpc预测控制算法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This program implements a controller that uses a planning
% strategy for the surge tank example.
%
% Kevin Passino
% Version: 4/19/01
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize variables

clear

% Set the length of the simulation

Nnc=300;

T=0.1;      % Sampling rate

% As a reference input we use a square wave （define one extra 
% point since at the last time we need the reference value at
% the last time plus one）

timeref=1:Nnc;
r（timeref）=3.25-3*square（（2*pi/150）*timeref）; % A square wave input
%r（timeref）=3.25*ones（1Nnc+1）-3*（2*rand（1Nnc+1）-ones（1Nnc+1））; % A noise input
ref=r（1:Nnc）;  % Then use this one for plotting

time=1:Nnc;
time=T*time; % Next make the vector real time


% We assume that the parameters of the surge tank （truth model） are:

abar=0.01;  % Parameter characterizing tank shape （nominal value is 0.01）
%abar=0.05;  % Parameter characterizing tank shape 
bbar=0.2;   % Parameter characterizing tank shape （nominal value is 0.2）
cbar=1;     % Clogging factor representing dirty filter in pump （nominally 1）
%cbar=0.8;     % Clogging factor representing dirty filter in pump （dirty case）
dbar=1;     % Related to diameter of output pipe 
g=9.8;      % Gravity

% Controller parameters

% Model used in planning strategy:

abarm=0.002;  % Parameter characterizing tank shape 
bbarm=0.2;   % Parameter characterizing tank shape 
cbarm=0.9;  % Clogging factor representing dirty filter in pump 
dbarm=0.8;   % Related to diameter of output pipe 

% Next use a set of preset controllers （think of each as a plan template）
% In this case we use PI controllers with a grid of Kp Ki gains around the 
% ones that we had designed for the PI controller for the **model**

% Guess at values 
Kp=0.01;
Ki=0.3;

%Kpvec=Kp; % Vector of gains in the region of Kp Ki
%Kivec=Ki;

Kpvec=0:0.05:0.2; % Vector of gains in the region of Kp Ki
%Kivec=Ki;
Kivec=0.15:0.05:0.4;

%Kpvec=0.2:0.2:1.8; % Vector of gains in the region of Kp Ki
%Kivec=0.05:0.01:0.11;

% Set weights of cost function used to select plans
w1=1;
w2=1;

% Set the length of time that will project into the future N
%NN=[20]; % To test just one value for N
NN=[151015172025303335363738394045505152535556596061656970];  % To test multiple values of the projection length
N=max（NN）; % For initialization and allocation of variables below


% First compare the truth model characteristics to the model

height=0:.1:10;

for i=1:length（height）
crosssect_truth（i1）=abs（abar*height（i）+bbar）;
crosssect_model（i1）=abarm*height（i）^2 +bbarm;
end

figure（1）
clf
plot（heightcrosssect_truth‘k-‘heightcrosssect_model‘k--‘）
grid
ylabel（‘Cross-sectional area‘）
xlabel（‘Height‘）
title（‘Cross-sectional area for plant （solid） and model （dashed）‘）


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next st