
% PROGRAMS used for paper tuningshort.tex:
% Written by S. Skogestad 02 Dec.01


%0. Data and tunings
     w= logspace(-4,4,400);
      numg=1;theta=0;deng=conv([1 2 1],[1 2 1]);

%1A. Find ZN-tunings
     Kc =1; taui=9999; taud=0; 
     num=conv(numg,conv([1 1/taui],[taud 1])); 
     den=conv(deng,conv([1 0],[0.000001*taud 1])); 
    [A,vA,w180,wc,GM,PM,Ms,Mt]=minbode(num,den,Kc,theta,w);
     Ku=GM, Pu = 2*pi/w180

     Kc = 0.45*Ku, taui = Pu/1.2            % ZN-PI
     Kc = 0.313*Ku, taui=2.2*Pu             % Tyreus-Luyben
     Kc = 0.6*Ku,  taui = Pu/2,   taud=Pu/8 % ZN-PID (ideal)
     Kc = 0.3*Ku,  taui = Pu/4,   taud=Pu/4 % ZN-PID (cascade)

% 1B. SIMC-tuning of first-order plus delay process

      k=1; theta=1; tau1=999;  
      tauc=theta % selected tuning parameter
      Kc = (1/k)*(tau1/(tauc+theta)), 
      tauI=min(tau1,4*(tauc+theta))

    num=[tauI 1];den=conv([tau1 1],[tauI 0]);
    [A,vA,w180,wc,GM,PM,Ms,Mt]=minbode(num,den,Kc*k,theta,w);


%2. Performance analysis using tunepid (simulation time 300, ode15s):
      Kc=0.3;taui=1.5;taud=0;
      tol=.1;
     tstop=300;
      td=400;tys=0.1;sim('tunepid');perf(1)=iae(Tid,y-ys);perf(2)=tv(u);
      td=0.1;tys=400;sim('tunepid');perf(3)=iae(Tid,y-ys);perf(4)=tv(u);perf

% 3.Stability analysis:
      [Gm,Pm,Ms,Mt]=marginsminbode(numg,deng,theta,Kc,taui,taud);
      [Gm,Pm,Ms,Mt]=marginss(numg,deng,Kc,taui,taud);         % NO DELAY: check with mu-toolbox

% 3a. combine stability and performance , Cascade PID:

      [Gm,Pm,Ms,Mt]=marginsminbode(numg,deng,theta,Kc,taui,taud);
      tstop=300;tol=0.1;td=400;tys=0.1;sim('tunepid');perf(1)=iae(Tid,y-ys);perf(2)=tv(u);
      td=0.1;tys=400;sim('tunepid');perf(3)=iae(Tid,y-ys);perf(4)=tv(u);
      margins=[Gm,Pm,Ms,Mt], perf

%3b. Ideal PID: Performance analysis using tunepid (simulation time 300, ode15s):
      corr = 1 + taud/taui; Kc=Kc*corr, taui=taui*corr, taud=taud/corr, % ZN cascade -> ideal

      [Gm,Pm,Ms,Mt]=marginsminbodeideal(numg,deng,theta,Kc,taui,taud);
      tstop=300;tol=0.1;td=400;tys=0.1;sim('tunepidideal');perf(1)=iae(Tid,y-ys);perf(2)=tv(u);
      td=0.1;tys=400;sim('tunepidideal');perf(3)=iae(Tid,y-ys);perf(4)=tv(u);
      margins=[Gm,Pm,Ms,Mt], perf

% 4. IAE-Optimal tuning for input (load) disturbance:
      global Kc taui taud;  
      td=0.1;tys=400; taud=0; xpi0= [Kc taui];  tol=.1;      xpi =  fminsearch(@funiaepi,xpi0) 
      td=0.1;tys=400; xpid0=[Kc taui taud];  tol=.1; xpid = fminsearch(@funiaepid,xpid0) 
% Med taui=taud
      td=0.1;tys=400; xpi0= [Kc taui];  tol=.1;      xpi =  fminsearch(@funiae,xpi0) 

% 4B. IAE-Optimal tuning with KI as a parameter (thought it would be better behaved, but does not seem to be the case):
      global Kc taui taud;  
      td=0.1;tys=400; taud=0; xpi0= [Kc/taui taui];  tol=.1;      xpi =  fminsearch(@funiaepiki,xpi0) 
      td=0.1;tys=400; xpid0=[Kc/taui taui taud];  tol=.1; xpid = fminsearch(@funiaepidki,xpid0) 
% Med taui=taud
      td=0.1;tys=400; xpi0= [Kc/taui taui];  tol=.1;      xpi =  fminsearch(@funiaeki,xpi0) 

% 4C. IAE-Optimal tuning with ideal PID:
      global Kc taui taud;  
      td=0.1;tys=400; xpid0=[Kc taui taud];  tol=.1; xpid = fminsearch(@funiaepidideal,xpid0) 


% 5. Constrained optimization to minimize IAE
% With inequality constraints on Ms at various frequencies + equality constraint on peak of Ms:
%  Conclusion: fmincomn is UNRELEIABLE - it even exceeds bounds in some cases! (negative Kc...?)
% PI
      global Kc taui taud numg deng theta;  td=0.1;tys=400;  tol=0.01; 
 A=[];b=[];Aeq=[];beq=[];lb=[0 0.1];ub=[10 999];x0=[Kc taui]; % pi
 x = fmincon(@funiaepiki,x0,A,b,Aeq,beq,lb,ub,@nonlconms)
% PID
  A=[];b=[];Aeq=[];beq=[];lb=[0 1 0];ub=[2 8 4];x0=[Kc taui taud];
  x = fmincon(@funiaepidki,x0,A,b,Aeq,beq,lb,ub,@nonlconms)
% PID with taui=taud
 A=[];b=[];Aeq=[];beq=[];lb=[0 1];ub=[2 8];x0=[Kc taui]; % pid
 x = fmincon(@funiaeki,x0,A,b,Aeq,beq,lb,ub,@nonlconms)

%6. Constrained optimization to minimize KI
      global Kc taui taud numg deng theta;  ki=Kc/taui; td=0.1;tys=400; tol=10; 
 A=[];b=[];Aeq=[];beq=[];lb=[0.0001 0.0001];ub=[100 100];x0=[Kc ki]; % pi
 xoptki = fmincon(@funkipi,x0,A,b,Aeq,beq,lb,ub,@nonlconmski)

% 7. We can also do other things (not in the paper): e.g. ISE-optimal PI for setpoints
      global Kc taui taud;  
      td=400;tys=0.1; xpi0= [Kc taui];  taud=0; tstop=300; tol=.1;  xpi =  fminsearch(@funisepi,xpi0) 

