Can anyone help me to rectify the error?

6 次查看（过去 30 天）
Sudipa
Sudipa 2021-1-19
回答： Anudeep Kumar 2025-3-14
Ran in:
main.m
clear
epsilon = 1/12;
gamma = 1/15;
beta = 0.2;
theta1 = 0.0082;
theta2 = 0.0289;
mu = 1.01/(70*365);
Lambda = 0.379/(70*365);
rho = 205;
delta0 = 1050;
r= 1/7;
r0 = 0.04;
delta1 = 0.0399;
d = 1/10;
psi = 0.0014;
zeta1 = 0;
zeta2 = 0.67;
alpha = 0.0238;
phi = 3;
A1 = 1;
A2 = 1;
B1 = 50;
B2 = 50;
T = 365; %final time
N = 200;
x0 = [100; 20; 20; 10; 1000; 20; 30]; %initial value of x
u0 = [0; 0]; %initial value of u
lambdafinal =[-0.1;0;0;-0.1;0;0;0]; %final values of lambda
t = linspace(0,T,N+1);
h = T/N;
x = zeros(7,N+1); %initialize x
u = zeros(2,N+1); %initialize u
lambda = zeros(7,N+1); %initialize x
x (:,1) = x0; %assign the x0
u (:,1) = u0; %assign the u0
lambda(:,N+1) = lambdafinal; %assign the lambdafinal
%R0
%1/2*(zeta1 * Lambda*d/(k1*k2) + zeta2 * delta1*epsilon/(gamma*k4))...
%+sqrt((1/2 *(zeta1*Lambda*d/(k1*k3)?zeta2*delta1*epsilon/(gamma*k4)))^2...
%+phi^2 * (beta*epsilon*(theta1*d+theta2*k3)*(delta0/(gamma?delta1)))/(gamma*
%k1*k3*k4*(rho/(mu-Lambda))))
parameters=[epsilon,gamma,beta,theta1,theta2,mu,Lambda,rho,delta0,r, r0, delta1, d, psi , zeta1, zeta2,alpha,phi,A1,A2,B1,B2];
k = 1;%counter of the iteration
delta = 0.001; %error bound
test = -1; %error
while ( test<0 && k<1000)
oldx = x;
oldu = u;
oldlambda = lambda;
%forward Runge?Kutta 4th with 3 input algorithm for state
for i = 1:N
k1(1:7,1) = state(t(i), x(:,i) ,u(:,i) ,parameters);
k2(1:7,1) = state(t(i)+h/2,x(:,i )+h*k1/2,(u(:,i)+u(:,i+1))/2,parameters);
k3(1:7,1) = state(t(i)+h/2,x(:,i )+h*k2/2,(u(:,i)+u(:,i+1))/2,parameters);
k4(1:7,1) = state(t(i)+h,x(:, i )+h*k3,u(:,i+1),parameters);
x (:,i+1) = x(:,i)+ (h/6)*(k1+2*k2+2*k3+k4);
end
%backward Runge?Kutta 4th with 3 input algorithm for adjoint
for i = 1:N
j = N+2-i;
k1(1:7,1) = adjoint(t(j),lambda(:,j), x(:,j),u(:,j),parameters);
k2(1:7,1) = adjoint(t(j)-h/2,lambda(:,j)-h*k1/2,(x(:,j)+x(:,j-1))/2 ,(u (:,j)+u(:,j-1))/2,parameters);
k3(1:7,1) = adjoint(t(j)-h/2,lambda(:,j)-h*k2/2,(x(:,j)+x(:,j-1))/2 ,(u (:,j)+u(:,j-1))/2,parameters);
k4(1:7,1) = adjoint(t(j)-h/2,lambda(:,j)-h*k3/2,x(:,j-1),u(:,j-1),parameters);
lambda(:,j-1) = lambda(:,j) - (h/6)*(k1+2*k2+2*k3+k4);
end
%Find controls
for i = 1:N+1
u(1, i ) = max(0,min(1,(beta*lambda(1,i)*(-x(1,i))*x(7,i)*phi +beta*lambda(2,i)*x(1,i)*x(7,i)*phi ...
+(lambda(6,i)-lambda(5,i))*x(5,i)*phi*(theta2*x(2,i) +theta1*x(3,i)))/(2*B1*(x(1,i)+x(2,i)+x(3,i)+x(4,i))))) ;
u(2, i ) = max(0,min(1,(lambda(3,i)-lambda(4,i))*r0*x(3,i)/(2*B2)));
end
%updates Control
c = 0.8;
u = (1-c)*u+c*oldu;
% for i=1:N+1
% if (u(1, i ) > oldu(1,i))
% u(1, i ) = (1 ? c) + oldu(1,i)*c;
% else
% u(1, i ) = oldu(1,i)*c;
% end
% if (u(2, i ) > oldu(2,i))
% u(2, i ) = (1 ? c) + oldu(2,i)*c;
% else
% u(2, i ) = oldu(2,i)*c;
% end
% end
%error bewteen old and new
tempu = min(delta*sum(abs(u),2)-sum(abs(oldu-u),2));
tempx = min(delta*sum(abs(x),2)-sum(abs(oldx-x),2));
templambda = min(delta*sum(abs(lambda),2)-sum(abs(oldlambda-lambda),2));
test = min(tempu,min(tempx,templambda));
%The cost
%trapz(A1*oldx(2,:)+A2*oldx(3,:)+B1*oldu(1,:).^2+B2*oldu(2,:).^2)...
% ?0.1*oldx(1,N+1)? 0.1*oldx(4,N+1)
%trapz(A1*x(2,:)+A2*x(3,:)+B1*u(1,:).^2+B2*u(2,:).^2) ...
% ?0.1*x(1,N+1)?0.1*x(4,N+1)
k=k+1;
end
plot(t,x(3,:),t,x(7,:))
title (' infectious ')
legend('host','vector')
state.m
function dxdt = state(t,x,u,parameters)
parameters = num2cell(parameters);
[epsilon,gamma,beta,theta1,theta2,mu,Lambda,rho,delta0,r, r0, delta1, d, psi, zeta1, zeta2,alpha,phi] = deal(parameters{:};
Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for mismatched delimiters.
dxdt = zeros(7,1);
dxdt(1) = rho + Lambda*(x(1)+x(2)+x(4)) + Lambda*(1-zeta1)*x(3)+psi*x(4) -beta* phi*x(7)*x(1)*(1-u(1))/(x(1)+x(2)+x(3)+x(4))-mu*x(1);
dxdt(2) = beta*phi*x(7)*x(1)*(1-u(1))/(x(1)+x(2)+x(3)+x(4))+ Lambda*zeta1*x(3)-d*x(2)-mu*x(2);
dxdt(3) = d*x(2) - (r + alpha + mu+r0*u(2))*x(3);
dxdt(4) = (r+r0*u(2))* x(3)- (psi + mu)*x(4);
dxdt(5) = delta0 + delta1*(x(5)+x(6)) + delta1*(1-zeta2)*x(7) - phi*theta1*x(3)*x(5)*(1-u(1))/(x(1)+x(2)+x(3)+x(4)) ...
- phi*theta2*x(2)*x(5)*(1-u(1))/(x(1)+x(2)+x(3)+x(4)) - gamma*x(5);
dxdt(6) = phi*theta1*x(3)*x(5)*(1-u(1))/(x(1)+x(2)+x(3)+x(4)) + phi*theta2*x(2)*x(5)*(1-u(1))/(x(1)+x(2)+x(3)+x(4))...
+ delta1*zeta2*x(7)-epsilon*x(6)-gamma*x(6);
dxdt(7) = epsilon*x(6)-gamma*x(7);
adjoint.m
function dlambdadt = adjoint(t,lambda,x,u,parameters)
parameters = num2cell(parameters);
[epsilon,gamma,beta,theta1,theta2,mu,Lambda,r, r0, delta1, d, psi, zeta1, zeta2,alpha,phi,A1,A2] = deal(parameters{:});
dlambdadt = zeros(7,1);
dlambdadt(1) = -(lambda(1)*(Lambda-mu+beta*(1-u(1))*x(1)*x(7)* phi/(x(1)+x(2)+x(3)+x(4))^2-beta*(1-u(1))*x(7)*phi/(x(1)+x(2)+x(3)+x(4))) +lambda(2)*(beta*(1-u(1))*x(7)*phi/(x(1)+x(2)+x(3)+x(4))-beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2 +lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2+theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2) +lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 -theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2);
dlambdadt(2) = -(A1+d*lambda(3) +lambda(2)*(-d-mu-beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2) +lambda(1)*(Lambda+beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2) +lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2+theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2-theta2*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4))) +lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 -theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2+theta2*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4))));
dlambdadt(3) = -(A2+lambda(3)*(-alpha-mu-r0*u(2)-r)+lambda(4)*(r0*u(2)+r) +lambda(1)*((1-zeta1)*Lambda+beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2) +lambda(2)*(zeta1*Lambda-beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2) +lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 +theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 -theta1*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4))) +lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+ x(4))^2 -theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2+theta1*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4)));
dlambdadt(4) = -(lambda(4)*(-mu-psi)+ lambda(1)*(Lambda+beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2+psi) -lambda(2)*(beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2) +lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 +theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2) +lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2) -theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2);
dlambdadt(5) = -(lambda(5)*(-gamma+delta1-theta2*(1-u(1))*x(2)*phi/(x(1)+x(2)+x(3)+x(4))-theta1*(1-u(1))*x(3)*phi/(x(1)+x(2)+x(3)+x(4))) +lambda(6)*(theta2*(1-u(1))*x(2)*phi/(x(1)+x(2)+x(3)+x(4))+theta1*(1-u(1))*x(3)*phi/(x(1)+x(2)+x(3)+x(4))));
dlambdadt(6) = -(lambda(6)*(-gamma-epsilon)+delta1*lambda(5)+lambda(7)*epsilon);
dlambdadt(7) = -(-gamma*lambda(7) +delta1*(1-zeta2)*lambda(5)+ delta1*zeta2*lambda(6) -beta*lambda(1)*(1-u(1))*x(1)*phi/(x(1)+x(2)+x(3)+x(4)) +beta*lambda(2)*(1-u(1))*x(1)*phi/(x(1)+x(2)+x(3)+x(4)));

请先登录，再回答此问题。

回答（1 个）

Anudeep Kumar
Anudeep Kumar 2025-3-14
Hey Sudipa,
I ran your code on my end and was able to reproduce the errors you encountered. Based on the R2024b definition of the functions you used; I identified two main issues:
Firstly, for line 3 in “state.m” -and “adjoint.m” the use of the MATLAB function “deal()” requires the number of input arguments to be equal to the number of output arguments. You were passing all 22 arguments in your variable “parameters” which was leading to an error.
Assuming you want to pass the first 18 arguments the following change in state.m:
>>[epsilon,gamma,beta,theta1,theta2,mu,Lambda,rho,delta0,r, r0, delta1, d, psi, zeta1, zeta2,alpha,phi] = deal(parameters{1:18});
And the below in adjoint.m:
>>[epsilon,gamma,beta,theta1,theta2,mu,Lambda,r, r0, delta1, d, psi, zeta1, zeta2,alpha,phi,A1,A2] = deal(parameters{1:18});
This should resolve the error.
If you need to pass specific arguments, you can use cell indexing or create a new variable accordingly.
Here's a link to guide you on accessing data in a cell array:
https://www.mathworks.com/help/matlab/matlab_prog/access-data-in-a-cell-array.html
Also, here’s the documentation for the deal() function:
https://www.mathworks.com/help/matlab/ref/deal.html
Secondly, in the definition of “adjoint.m”, there was an imbalance in the number of opening and closing parentheses. Assuming the formula was missing parentheses at the end, the following changes resolved the issue:
At line 6:
>> dlambdadt(1) = -(lambda(1)*(Lambda-mu+beta*(1-u(1))*x(1)*x(7)* ...
phi/(x(1)+x(2)+x(3)+x(4))^2-beta*(1-u(1))*x(7)*phi/(x(1)+x(2)+x(3)+x(4)))...
+lambda(2)*(beta*(1-u(1))*x(7)*phi/(x(1)+x(2)+x(3)+x(4))...
-beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2 ...
+lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 ...
+theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2) ...
+lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 ...
-theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2)));
And at line 8 :
>>dlambdadt(3) = -(A2+lambda(3)*(-alpha-mu-r0*u(2)-r)+lambda(4)*(r0*u(2)+r)...
+lambda(1)*((1-zeta1)*Lambda+beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2)...
+lambda(2)*(zeta1*Lambda-beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))^2)...
+lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 ...
+theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 ...
-theta1*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4))) +lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+ x(4))^2 ...
-theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))^2 ...
+theta1*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4))));
If the formula is different, please review and ensure that all parentheses are correctly balanced.
I hope this helped resolve your errors!

类别

在 Help Center 和 File Exchange 中查找有关 MATLAB 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by