optfunf.m

Software for design and tuninig of SISO and MIMO contol systems

function [results, refreq]=OptFunf(Epsilon,fun,decade,numpdec,vub,vlb,x,y,MoM,p,m,q,qf,...
   qd,pd,Mp,Stop,Acc,Tcanc,order,mq);

% OPTFUN finds the maximum or minumum of the complementary sensitivity function
% subject to parametric process uncertanty over a range of frequencies
% and a given filter time constant (tuning parameter).

%       [results, refreq]=OptFun(Epsilon,fun,decade,numpdec,vub,vlb,x,MoM,p,m,q,Skip)
%       starts at x, initial guess of the worst-case plant parameters (x = [p1 p2 p3
%       p4 .... pn]), and finds the constrained optimum to the complementary
%       sensitivity function which is described in 'fun' (usually an
%       M-file: fun.m). The optimization is done over a frequency range
%       specified by the lower and upper decade: decade=[UpperD LowerD].
%       The lower and upper bounds of the process parameters are specified
%       using two vectors: vlb=[p1- p2- ... pn-] and vub=[p1+ p2+ ... pn+].
%       An initial value of the filter-time constant, Epsilon, must be specified.
%       The M-file fun.m should receive the parameters x, w, and Epsilon, and
%       should return a scalar value of the function to be maximized.

% INITIALIZE PARAMETERS
global xwclb wlb
clear results Optimum;
if isempty(x) | length(x)~=length(vub)
   x=(vub+vlb)/2;
end
% ------------------------- DEFINE PARAMETERS -------------------------------
npoints=round((decade(2)-decade(1))*numpdec);
refreq=[10^decade(1) 10^decade(2)];
freq=logspace(decade(1),decade(2),npoints);
results=zeros(npoints,2+length(x));

if strcmp(fun,'imc2sen') | strcmp(fun,'imc2com')| strcmp(fun,'imc2pseudo')
   [qd,qdnum,qdden]=qd_mat(Tcanc(2,:),mq,Epsilon(2),order,y);
end
options=optimset('LargeScale','off','Display','off');

% ---------------------------------------------------------------------------

% ----- COMPUTATION OF UPPER OR LOWER BOUND ON CLOSED-LOOP BODE DIAGRAM -----
fprintf('\n       Iteration             Magnitude         Frequency \n');
fprintf('       ---------             ---------         --------- \n');
for i=1:npoints
   [x,temp]=fmincon(fun,x,[],[],[],[],vlb,vub,[],options,y,Epsilon,MoM,freq(i),0,0,p,m,q,qf,qd,pd);
   Optimum=MoM*temp;
   results(i,:)=[freq(i) Optimum x];
   if MoM==-1
   	fprintf('          %3.0f                %.4f            %g\n',i,Optimum,freq(i));
   end
   if Stop>0 & Optimum<Stop
     break;	     
   end
end
results=results(1:i,:);lb=0;
[row,col]=size(results);
if MoM==1
	% Globality Check by backward re-computing of the lower bound //by Tom Feb 3,00
	for i=row:-1:1 % backward
      if i > 1 & any(results(i-1,3:end)~=results(i,3:end))
         temp = feval(fun,results(i,3:end),y,Epsilon,1,freq(i-1),0,0,p,m,q,qf,qd,pd);
         if temp < results(i-1,2)
            results(i-1,:)=[freq(i-1) temp results(i,3:end)];
         end
      end
      if abs(results(i,2)-0.7) < abs(lb-0.7)
         xwclb=results(i,3:end);
         wlb=results(i,1);
         lb=results(i,2);
      end
      
      fprintf('          %3.0f                %.4f            %g\n',...
      	row-i+1,results(i,2),freq(i));
	end
 	lwb=results;
%  	save graph4 fun Epsilon refreq lwb
else 
	% ---  Trying to insert the real optima of the upper bound
	% --------------------
	nmax=0;
	locmax=[];
	for i=2:row-1
 		if results(i,2)>results(i-1,2) & results(i,2)>results(i+1,2)
  			nmax=nmax+1;
  			locmax(nmax,:)=[i results(i,:)];
 		end
	end
	Utemp=vub; Ltemp=vlb;
	if nmax > 0
 		for i=1:nmax
  			n=col-1;  Xo=locmax(i,4:col+1);  Xo(n)=locmax(i,2);
  			Ltemp(n)=results(locmax(i,1)-1,1); Utemp(n)=results(locmax(i,1)+1,1);
        	[Xo,temp]=fmincon(fun,Xo,[],[],[],[],Ltemp,Utemp,[],options,y,Epsilon,-1,0,n,...
           Mp,p,m,q,qf,qd,pd);
  			Optimum=MoM*temp;
  			if Optimum > locmax(i,3)
   			results(locmax(i,1),:)=[Xo(n) Optimum Xo(1:n-1)];
  			end
 		end  % for
	end  % nmax>0
 	upb=results;
end