%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%?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