Error using + Matrix dimensions must agree.

조회 수: 2 (최근 30일)
senadjki mohamed islam
senadjki mohamed islam 2020년 7월 3일
답변: Are Mjaavatten 2020년 7월 5일
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

답변 (1개)

Are Mjaavatten
Are Mjaavatten 2020년 7월 5일
In line 43, x(:,i) is a column vector, while k1' is a row vector. Summing a column vector and a row vector is not allowed in older versions of Matlab. In Matlab 2014b:
>> a = [1;2;3];b = [10,20,30];
>> a + b
Error using +
Matrix dimensions must agree.
Your code runs without error in version 2020a. Here vector a is added to each element in b, resulting in a 3 by 3 array:
>> a = [1;2;3];b = [10,20,30];
>> a+b
ans =
11 21 31
12 22 32
13 23 33
But this is proably not what you want. Your glucoz function expects a three element vector, so you should simply use the column vector k1 instead of the transpose k1'.

카테고리

Help CenterFile Exchange에서 Language Fundamentals에 대해 자세히 알아보기

제품


릴리스

R2015a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by