clear all 
clc
global data
load Stationary
load data
OPz=30;
index=find(abs(z-OPz)==min(abs(z-OPz)));
P1(index)
mG_in=(6.1717e-5)*1.5;
P_feed=0.31e5+1.013e5;
P0=1.013e5;
u=[z(index)*0.01;mG_in;P_feed;P0/1e5];
%Linearizing
[A,B,C,D,E,F]=linearize(@slug_3D_lin,[],xs(1:3,index),[],u(1),u(2:4));
G=ss(A,B,C,D);
G_tf=tf(G);

GP1=G_tf(1);
GP1=zpk(GP1);
GP2=G_tf(2);
GP2=zpk(GP2);
GpT=G_tf(3);
GpT=zpk(GpT);
GW=G_tf(4);
GW=zpk(GW);
GQ=G_tf(5);
GQ=zpk(GQ);

[num_GP1,den_GP1]=tfdata(GP1,'v');
[num_GP2,den_GP2]=tfdata(GP2,'v');
[num_GpT,den_GpT]=tfdata(GpT,'v');
[num_GW,den_GW]=tfdata(GW,'v');
[num_GQ,den_GQ]=tfdata(GQ,'v');
A=10;
B=10;
% Kc=-25;%bar^-1;
% taui=10;%s
% Gc=Kc*tf([taui 1],[taui 0]);
% S=inv(1+GP1*Gc);
% T=1-S;

%PI regulator m/prosess variabel P1


Kc=-randperm(400);
taui=randperm(800);
for i=1:100
K=Kc(i);
tau=taui(i);
Gc=K*tf([tau 1],[tau 0]);
S=inv(1+GW*Gc);
T=1-S;
MS=norm(S,inf);
MT=norm(T,inf);
if ((MS+MT)<(A+B)) & (MS>1 & MT>1)
    A=MS;
    B=MT;
    a=K;
    b=tau;
else
end

end
SUM_M=A+B,MS=A,MT=B,K=a,tau=b




Kc=-randperm(200);%bar^-1;
taui=randperm(800);%s
Gc=Kc*tf([taui 1],[taui 0]);
S=inv(1+GP1*Gc);
T=1-S;
MS=norm(S,inf)
MT=norm(T,inf)
sum=MS+MT