msf_adj.m

Software for design and tuninig of SISO and MIMO contol systems

function [e]=msf_adj(pnummt,pdenmt,pdeadmt,mnummt,mdenmt,mdeadmt,pdnummt,pddenmt,cnummt,cdenmt,...
   order,x,y,Epsilon,subs,slbs,d_through_p,cnums,cdens,e_step)
%Model
%-----------------------------------------------
[Ksp,K]=msfgain(cnums,cdens,Epsilon(1));
[numMs]=mt2poly(mnummt,y);
[denMs]=mt2poly(mdenmt,y);
[mdeadx,mnumT,mdenT]=mt2dt(mdeadmt,y);

%Process
%-----------------------------------------------
if isempty(pnummt) | isempty(pdenmt)
   numPs=numMs;
   denPs=denMs;
   pdeadx=mdeadx;
   pnumT=mnumT;
   pdenT=mdenT;
else
   [numPs]=mt2poly(pnummt,x);
	[denPs]=mt2poly(pdenmt,x);
	[pdeadx,pnumT,pdenT]=mt2dt(pdeadmt,x);
end

denF=[];
for i=1:order
 denF(i,1)=Epsilon(1);
 denF(i,2)=1;
end
[r,c]=size(denF);
numFs=denF(1,:);
for i=2:r
 numFs=conv(numFs,denF(i,:));
end
if Epsilon==0
 numFs=1;
end

denF=[];
for i=1:order
 denF(i,1)=Epsilon(1);
 denF(i,2)=1;
end
[r,c]=size(denF);
denFs=denF(1,:);
for i=2:r
 denFs=conv(denFs,denF(i,:));
end
if Epsilon==0
 denFs=1;
end

%IMC Controller
denCs=mt2poly(cnummt,y);
numCs=mt2poly(cdenmt,y);

%Approximated feed-back form controller
%[Csnum,Csden]=pid_tf(numMs,denMs,mnumT,mdenT,denCs,numCs,denFs);

denCs=conv(denCs,denFs);


% Disturbance lag------------------
numPd=mt2poly(pdnummt,x);
denPd=mt2poly(pddenmt,x);


%Delays and Saturation bounds
%---------------------------------------------------------

Mdelay=mdeadx;
Pdelay=pdeadx;
Satub=str2num(subs);
Satlb=str2num(slbs);
dimn=length(numMs);
dimd=length(denMs);
numMs=numMs/denMs(dimd);
denMs=denMs/denMs(dimd);
Md=denMs;
Md(dimd+1:9)=0;
Mnt=[];
for i=1:(dimd-1)
	if (dimn-i+1)>0;
		Mnt(dimd-i)=numMs(dimn-i+1);
	else (dimn-i+1)<=0;
		Mnt(dimd-i)=0;
	end
end
Mnt(dimd:8)=0;
if length(K)<dimd;
   Ktemp=K;
   for i=1:dimd;
      if i <=length(K);
         K(1,dimd+1-i)=Ktemp(1,length(K)+1-i);
      else
         K(1,dimd+1-i)=0;
      end
   end
end
K(dimd+1:9)=0;  

%load msfgain.mat
if Pdelay==0
	Pdelay=0.0001;
end
if Mdelay==0
	Mdelay=0.0001;
end
time=50;
TimeScale=200;
%Mnum=Xmnums{1,1};
%Mden=Xmdens{1,1};
e=Epsilon(1);
f=crfig(3,5,634,410,'MSF_Filter.T.C Modifying Process','k','figure','on');
whitebg('white');
%if length(x)~=0;
%	store=x;
%end;
%x=(vub+vlb)/2;
%e_step=0.1;
options = simset('SrcWorkspace','current');
Setpoint=1;
Disturb=0;
[t]=sim('msf1df_unsat',time,options);
ceffort_min=min(simout(:,2));
output_max=max(simout(:,1));
%e_step=0.05;

while ceffort_min-Satlb < 0;
	e=e+e_step;
	Ksp=Md(1)/e^(dimd-dimn)/numMs(1);
	filter=1;
	for i=1:dimd-dimn
		filter=conv(filter,[e 1]);
   end
   K=Ksp*(conv(numMs,filter))-Md(1:dimd);
	K(dimd+1:9)=0;  
	[t]=sim('msf1df_unsat',time,options);
	ceffort_min=min(simout(:,2));
	output_max=max(simout(:,1));
	if ceffort_min-Satlb > 0; 
		break,
	end
end
%	save msfgain_adj Ksp K e x;
	axis('auto');
	subplot(2,1,1),plot(t,simout(:,1),'m-');	
	c=axis;
	axis([c(1) time c(3) c(4)]);
	title(['MSF_Saturation control effort when Epsilon = ',num2str(e)])
	xlabel(' '); ylabel('Output'); grid on;
	hold on;
	subplot(2,1,2),plot(t,simout(:,2),'m-',[min(t) max(t)],[Satlb Satlb],'m-.');	
	c=axis;
	axis([c(1) time c(3) c(4)]);
	title(['MSF_Saturation control effort when Epsilon = ',num2str(e)])
	xlabel('Time '); ylabel('Control effort'); grid on;
	hold on;