close all
clear all
clc
G0=300; %initial positions
I0=2.5;
q=1/1;
A=90;
D(1)=A;
Gd=80;
Z(1)=0;
Ib=2.5;
Comm=0;
x(:,1)=[G0 0 20]'; %states initial values
xob(:,1)=[G0 0 20]';
lam(:,1)=[0 0 0]';
% patient 1
p1=0;
p2=0.0142;
p3=154*10^-5;
n=0.2814;
Pn=[p1 p2 p3 n];
%%% patient 2
%p1=0;
%p2=0.0072;
%p3=2.16*10^-5;
%n=0.2865;
Pn=[p1 p2 p3 n];
Comande(1)=10;
i=1;
B=0.05;
t(1)=0;
h=1/6;%0.005;1/6
% ODE solving
% opt1=odeset('RelTol',1e-10,'AbsTol',1e-20,'NormControl','off');
%[T,X] = ode45(@(t,x) glucoz(t,x,Pn,Gb,Comm),Ts,x0);
%V(1)=0;
YY3(1)=0;
DD(1)=0;
a=1;
while t<18000
D(i+1)=A*exp(-B*i*h);
k1=h*glucoz(x(:,i),Pn,Gd,Comande(i),D(i));
k2=h*glucoz(x(:,i)+k1'/2,Pn,Gd,Comande(i),D(i)); %the error in the + signe
k3=h*glucoz(x(:,i)+k2'/2,Pn,Gd,Comande(i),D(i));
k4=h*glucoz(x(:,i)+k3',Pn,Gd,Comande(i),D(i));
x(:,i+1)=x(:,i)+(1/6)*(k1+2*k2+2*k3+k4);
DD(i+1)=-A*B*exp(-B*i*h);
D2D(i+1)=A*B^2*exp(-B*i*h);
Food=0;
if i*h>Food
D(i+1)=2*A*exp(-B*(i*h-Food));
DD(i+1)=-A*B*exp(-B*(i*h-Food));
D2D(i+1)=A*B^2*exp(-B*(i*h-Food));
else
D(i+1)=D(i+1);
DD(i+1)=DD(i+1);
D2D(i+1)=D2D(i+1);
end
%%% tracking error %%%%%%%%%
e= x(1,i+1)-Gd;%0.1*sign(x(3,i+1)-Ib);
%%%%%%%%%%%%%%% First derivative %%%%%
edot=-x(1,i+1)*x(2,i+1)+D(i);
%%%%%%%%%%%% second derivative %%%%%%%%%%%
e2dot=-x(1,i+1)*x(2,i+1)^2+D(i+1)*x(2,i+1)+p2*x(1,i+1)*x(2,i+1)-p3*(x(3,i+1)-Ib)*x(1,i+1)+DD(i+1);
%%%%%%%%%%%%%%% Third Derivative %%%%%%%%%%%%%
%f2=[D(i+1)-x(1,i+1)*x(2,i+1)]*x(2,i+1)^2-2*p3*x(1,i+1)*x(2,i+1)*x(3,i+1)-3*p3*Ib*x(1,i+1)*x(2,i+1)-3*p2*x(1,i+1)*x(2,i+1)^2+x(2,i+1)*DD(i+1)+p3*x(3,i+1)*D(i+1)+p2*p3*x(1,i+1)*x(3,i+1)-p2*p3*Ib*x(1,i+1)-p2^2*x(1,i+1)*x(2,i+1)+D2D(i+1)+p3*x(3,i+1)*D(i+1)+p3*n*(x(3,i+1)-Ib)*x(2,i+1);
F=[D(i+1)-x(1,i+1)*x(2,i+1)]*[p2*x(2,i+1)-p3*(x(3,i+1)-Ib)+x(2,i+1)^2]-[p3*(x(3,i+1)-Ib)-p2*x(2,i+1)]*(D(i+1)+p2*x(1,i+1)-2*x(1,i+1)*x(2,i+1))+DD(i+1)*x(2,i+1)+D2D(i+1)+p3*n*x(1,i+1)*(x(3,i+1)-Ib);
%e3dot=F-p3*x(1,i+1)*U;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% sliding surface %%%%%%%%%%
lambda1=0.01;%0.01;%0.01
lambda2=0.0001;%0.0001;%0.0001
QQ=.000;
KK=0.5;%0.5; 0.1 with Q=0.001, Patien 1
S=(q^2*e2dot+q*lambda1*edot+lambda2*e);%+0.001*sign(x(3,i+1)-Ib);KK=0.5;QQ=0.1;
surface(i)=S;
if KK>0
KK=KK+h*(abs(S)*sign(abs(S)-0.0001));
else
KK=KK;
end
%KK=0.9;
%KK=abs(p3*x(1,i+1)*10*F/(F^2+0.1));
YY=1.7159*[exp(a*S)-1]/[exp(a*S)+1];
g(i)=p3*x(1,i+1)*q^3;
GG=x(1,i+1)*p3;
U=((q^3)*F+(q^2)*lambda1*e2dot+q*lambda2*edot+KK*YY+QQ*S)/GG;
%___________________________________________________________
Comande(i+1)=U;
X=x(2,i+1);
I=x(3,i+1);
t(i+1)=i*h;
temps(i)=i*h;
i=i+1;
end
% Output
%torque inputs computation from the 7th,8th states inside ODE
%tt=0:(T(end)/(length(Comm)-1)):T(end);
%theta1 error plot
figure(1)
plot(t/60,x(1,:),'k')%,t/60,xob(1,:))
grid
title('G')
ylabel('error ')
xlabel('time (min)')
figure(2)
plot(t/60,Comande),
title('Commande')
figure(3)
plot(t/60,x(3,:),'k')
grid
title('X3')
ylabel('Insuline in plasma')
xlabel('time (min)')
figure(4)
plot(t/60,x(2,:),'k')
grid
title('X2')
ylabel('insuline effect')
xlabel('time (min)')
%theta2 error plot
figure(5)
plot(temps/60,surface)
grid
title('surfaces')
ylabel('insuline effect')
xlabel('time (min)')
%theta2 error plot
