clear
h=0.1;t=0.1;N=1500;
r1=0.1;r2=4;r3=4;r4=4;r5=4;r6=4;
yr0=0;e10=0;Yp0=0;
a=4;b=5;c=6;km=40;d=3;
am0=a*b*c;am1=a*b+a*c+b*c;am2=a+b+c;bm0=km*d;bm1=km;
T1=1;T2=2;T3=4;k=0.4;
ap1=1/(T1*T2);ap2=1/T1+1/T2;b0=1/T3;kp=k*T3/(T1*T2);
a0=6;g0=d*a0;g1=d+a0;
xp10=0;xp20=0;xp30=0;xm10=0;xm20=0;xm30=0;
k00=3.5;c10=-0.5;c20=3.1;d00=-3.5;d10=0.05;d20=-0.63;
q10=0;q20=0;q30=0;q40=0;q50=0;q60=0;
v110=0;v120=0;v210=0;v220=0;u0=0;
%yr=[ones(1,250),zeros(1,250),ones(1,250),zeros(1,250),ones(1,250),zeros(1,250)];
yr=[ones(1,1500)];
for i=1:N
xm1=xm10+h*xm20;xm2=xm20+h*xm30;xm3=xm30+h*(-am0*xm10-am1*xm20-am2*xm30+yr0);
ym(i)=bm0*xm1+bm1*xm2;
xp1=xp10+h*xp20;xp2=xp20+h*xp30;xp3=xp30+h*(-ap1*xp20-ap2*xp30+u0);
Yp=kp*b0*xp1+kp*xp2;
yp(i)=Yp;
if i>=300
yp(i)=Yp+0.3;
else
yp(i)=Yp;
end
e1=yp(i)-ym(i);E1(i)=e1;
v11=v110+h*v120;v12=v120+h*(u0-g0*v110-g1*v120);
v21=v210+h*v220;v22=v220+h*(Yp0-g0*v210-g1*v220);
q1=q10+h*(-4*q10+yr0);
q2=q20+h*(-4*q20+v110);
q3=q30+h*(-4*q30+v120);
q4=q40+h*(-4*q40+Yp0);
q5=q50+h*(-4*q50+v210);
q6=q60+h*(-4*q60+v220);
Q1(i)=q1;Q2(i)=q2;Q3(i)=q3;Q4(i)=q4;Q5(i)=q5;Q6(i)=q6;
k0=k00+h*(-r1*q10*e10);
c1=c10+h*(-r2*q20*e10);
c2=c20+h*(-r3*q30*e10);
d0=d00+h*(-r4*q40*e10);
d1=d10+h*(-r5*q50*e10);
d2=d20+h*(-r6*q60*e10);
K0(i)=k0;C0(i)=c1;C1(i)=c2;D0(i)=d0;D1(i)=d1;D2(i)=d2;
u=k0*yr(1,i)+c1*v11+c2*v12+d0*Yp+d1*v21+d2*v22-(r1*q1^2+r2*q2^2+r3*q3^2+r4*q4^2+r5*q5^2+r6*q6^2)*e1;
U(i)=u;
xm10=xm1;xm20=xm2;xm30=xm3;
xp10=xp1;xp20=xp2;xp30=xp3;
k00=k0;c10=c1;c20=c2;d00=d0;d10=d1;d20=d2;
q10=q1;q20=q2;q30=q3;q40=q4;q50=q5;q60=q6;
e10=e1;Yp0=Yp;u0=u;yr0=yr(1,i);
v110=v11;v120=v12;v210=v21;v220=v22;
tt(i)=t;t=t+h;
end
figure(1);
plot(tt,yp,'-r')
h=legend('yp','yr',4);
xlabel('t(s)');
ylabel('yp');
figure(2);
plot(tt,E1);
xlabel('t(s)');
ylabel('E1');
figure(3);
plot(tt,U,'m');
xlabel('t(s)');
ylabel('U');