查看: 146|回复: 3|关注: 0

[已答复] matlab代码纠正错误

[复制链接]

新手

10 麦片

财富积分


050


1

主题

2

帖子

0

最佳答案
发表于 2020-7-1 10:45:42 | 显示全部楼层 |阅读模式
我写的代码如下:function ODEfunctiongroup_boke_wo
%  动力学ODE方程模型的参数估计
%  此例数据只有t,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10
%   k1-k12 为待拟合参数
%   The variables  y here are y(1)=x1, y(2)=x2, y(3)=x3,y(4)=x4,y(5)=x5,y(6)=x6,y(7)=x7,y(8)=x8,y(9)=x9,y(10)=x10.
%
clear all;
clc;
global t0
k0 = [5 5 5 5 5 5 5 5 5 5 5 5];        % 参数初值
lb = [0 0 0 0 0 0 0 0 0 0 0 0];        % 参数下限
ub = [+inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf];    % 参数上限
x0 = [2276.45 20.11 0.029 0.0010 0.19 0.23 0.12 0.045 119.24 0.089];  %初始状态 %动力学数据:
% t      x1    x2    x3   x4   x5      x6  x7   x8      x9   x10
ExpData = ...
[0 2276.45  20.11  0.03  0.00  0.19  0.23  0.12  0.04  119.24  0.09
4 2003.47  67.60  0.22  0.00  0.57  0.19  0.18  0.29  746.86  0.16
8 2032.79  77.54  0.29  0.00  0.93  0.15  0.24  0.50  1101.89 0.23
12 2149.09  56.71  0.21  0.00  0.94  0.15  0.24  0.47  1176.67 0.35
16 2011.55  42.71  0.15  0.01  0.89  0.14  0.26  0.47  995.61  0.33
20 1839.29  48.31  0.15  0.01  0.79  0.13  0.27  0.56  864.48  0.52
24 2003.07  34.93  0.19  0.02  0.88  0.11  0.28  0.79  737.48  0.60
28 1923.94  37.12  0.16  0.02  1.07  0.13  0.30  0.96  670.42  0.60
32 1790.84  60.29  0.19  0.03  1.34  0.14  0.52  1.54  672.09  1.16
36 1909.44  45.17  0.19  0.03  1.25  0.13  0.48  0.90  710.46  1.29
40 1700.47  45.07  0.21  0.03  1.06  0.14  0.33  0.44  672.98  0.78
];
t0 = ExpData(:,1);  
yexp = ExpData(:,2:11);                 % yexp: 对应实验数据[x1 x2 x3 x4 x5 x6 x7 x8 x9 x10]
%% 使用函数fmincon()进行参数估计
% [k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
% fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
% fprintf('\tk1 = %.4f\n',k(1))
% fprintf('\tk2 = %.4f\n',k(2))
% fprintf('\tk3 = %.4f\n',k(3))
% fprintf('\tk4 = %.4f\n',k(4))
% fprintf('\tk5 = %.4f\n',k(5))
% fprintf('\tk6 = %.4f\n',k(6))
% fprintf('\tk7 = %.4f\n',k(7))
% fprintf('\tk8 = %.4f\n',k(8))
% fprintf('\tk9 = %.4f\n',k(9))
% fprintf('\tk10 = %.4f\n',k(10))
% fprintf('\tk11 = %.4f\n',k(11))
% fprintf('\tk12 = %.4f\n',k(12))
% fprintf('The sum of the squares is: %.1e\n\n',fval)
% k_fmincon = k;
% 这一步通常被省略,通过反复迭代初始值得到最优解,加上后可以降低对初始值的依赖。
% 以函数fmincon()估计得到的结果为初值,使用函数lsqnonlin()进行参数估计
% 需要指出,这种方法并非在所有场合均有效,但有时确实可以改善求解效果。
%% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
   lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
Output(k,ci,resnorm)
%% 以函数fmincon()估计得到的结果为初值,使用函数lsqnonlin()进行参数估计
% k0 = k_fmincon;
% [k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
%    lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);      
% ci = nlparci(k,residual,jacobian);
% fprintf('\n\n以fmincon()的结果为初值,使用函数lsqnonlin()估计得到的参数值为:\n')
% Output(k,ci,resnorm)
%% ----------------------------------------------------
tspan = t0';     % 指定点时间,也可增加步数,只需多调用一次ode45函数
[t x] = ode45(@KineticEqs,tspan,x0,[],k);  
y(:,1:10) = x(:,1:10);
figure(1)
plot(tspan,y(:,1),'b',ExpData(:,1),yexp(:,1),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y1(x1)');
figure(2)
plot(tspan,y(:,2),'b',ExpData(:,1),yexp(:,2),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y2(x2)');
figure(3)
plot(tspan,y(:,3),'b',ExpData(:,1),yexp(:,3),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y3(x3)');
figure(4)
plot(tspan,y(:,4),'b',ExpData(:,1),yexp(:,4),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y4(x4)');
figure(5)
plot(tspan,y(:,5),'b',ExpData(:,1),yexp(:,5),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y5(x5)');
figure(6)
plot(tspan,y(:,6),'b',ExpData(:,1),yexp(:,6),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y6(x6)');
figure(7)
plot(tspan,y(:,7),'b',ExpData(:,1),yexp(:,7),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y7(x7)');
figure(8)
plot(tspan,y(:,8),'b',ExpData(:,1),yexp(:,8),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y8(x8)');
figure(9)
plot(tspan,y(:,9),'b',ExpData(:,1),yexp(:,9),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y9(x9)');
figure(4)
plot(tspan,y(:,10),'b',ExpData(:,1),yexp(:,10),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y10(x10)');
end
%% ------------------------------------------------------------------
function Output(k,ci,resnorm)
% Output.m
fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1))
fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2))
fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3))
fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4))
fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5))
fprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6))
fprintf('\tk7 = %.4f ± %.4f\n',k(7),ci(7,2)-k(7))
fprintf('\tk8 = %.4f ± %.4f\n',k(8),ci(8,2)-k(8))
fprintf('\tk9 = %.4f ± %.4f\n',k(9),ci(9,2)-k(9))
fprintf('\tk10 = %.4f ± %.4f\n',k(10),ci(10,2)-k(10))
fprintf('\tk11 = %.4f ± %.4f\n',k(11),ci(11,2)-k(11))
fprintf('\tk12 = %.4f ± %.4f\n',k(12),ci(12,2)-k(12))
fprintf('The sum of the squares is: %.1e\n\n',resnorm)
end
%% ------------------------------------------------------------------
% function f = ObjFunc4Fmincon(k,x0,yexp)
% global t0
% tspan = t0';
% [t x] = ode45(@KineticEqs,tspan,x0,[],k);  
% y(:,1:10) = x(:,1:10);      % 对应实验数据  x1 x2 x3 x4 x5 x6 x7 x8 x9 x10  
% f = sum((y(:,1)-yexp(:,1)).^2) + sum((y(:,2)-yexp(:,2)).^2)  ...
%     + sum((y(:,3)-yexp(:,3)).^2) + sum((y(:,4)-yexp(:,4)).^2)+sum((y(:,5)-yexp(:,5)).^2)+sum((y(:,6)-yexp(:,6)).^2)
%     +sum((y(:,7)-yexp(:,7)).^2)+sum((y(:,8)-yexp(:,8)).^2)+sum((y(:,9)-yexp(:,9)).^2)+sum((y(:,10)-yexp(:,10)).^2);   % 计算平方和,供fmincon调用
% end
%% ------------------------------------------------------------------
function f = ObjFunc4LNL(k,x0,yexp)
global t0
tspan = t0';
[t x] = ode45(@KineticEqs,tspan,x0,[],k);  
y(:,1:10) = x(:,1:10);            % 对应实验数据  x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
f1 = y(:,1) - yexp(:,1);
f2 = y(:,2) - yexp(:,2);
f3 = y(:,3) - yexp(:,3);
f4 = y(:,4) - yexp(:,4);
f5 = y(:,5) - yexp(:,5);
f6 = y(:,6) - yexp(:,6);
f7 = y(:,7) - yexp(:,7);
f8 = y(:,8) - yexp(:,8);
f9 = y(:,9) - yexp(:,9);
f10 = y(:,10) - yexp(:,10);
f = [f1; f2; f3; f4;f5;f6;f7;f8;f9;f10]; % 理论值与实验值的差值,残差
% figure(5)
% subplot(2,2,1)
% plot(f1);
% subplot(2,2,2)
% plot(f2);
% subplot(2,2,3)
% plot(f3);
% subplot(2,2,4)
% plot(f4);
end
%% ------------------------------------------------------------------
function dxdt = KineticEqs(t,x,k)   % 微分方程组
dxdt =  ...
[(-k(1)*x(1));
  (k(1)*x(1)-k(2)*x(3)-k(5)*x(3)-k(7)*x(3)-k(9)*x(3));
  (k(2)*x(3)-k(3)*x(4)+k(13)*x(5));
  (k(7)*x(3)-k(8)*x(5)-k(12)*x(5)-k(13)*x(5));
  (k(9)*x(3)-k(10)*x(6)-k(11)*x(6)-k(14)*x(6));
  (k(5)*x(3)-k(6)*x(7)+k(14)*x(6));
  (k(8)*x(5)+k(11)*x(6));
  (k(3)*x(4));
  (k(12)*x(5)+k(10)*x(6));
  (k(6)*x(7))
   ];
end
% code end

错误提示:索引超出数组元素的数目(12)。
               出错 ODEfunctiongroup_boke_wo>KineticEqs (line 159)
               (k(2)*x(3)-k(3)*x(4)+k(13)*x(5));
               出错 odearguments (line 90)
               f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.
               出错 ode45 (line 115)
               odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
               出错 ODEfunctiongroup_boke_wo>ObjFunc4LNL (line 131)
                [t x] = ode45(@KineticEqs,tspan,x0,[],k);
               出错 lsqnonlin (line 215)
               initVals.F = feval(funfcn{3},xCurrent,varargin{:});
                出错 ODEfunctiongroup_boke_wo (line 54)
                lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);
原因:  Failure in initial objective function evaluation. LSQNONLIN cannot continue.
回复主题 已获打赏: 0 积分

举报

论坛优秀回答者

6

主题

1440

帖子

306

最佳答案
  • 关注者: 70
发表于 2020-7-2 13:45:50 | 显示全部楼层
仅供参考,将
  1. k0 = [5 5 5 5 5 5 5 5 5 5 5 5];        % 参数初值
  2. lb = [0 0 0 0 0 0 0 0 0 0 0 0];        % 参数下限
  3. ub = [+inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf];    % 参数上限
复制代码

修改为
  1. k0 = [5 5 5 5 5 5 5 5 5 5 5 5 5 5];        % 参数初值
  2. lb = [0 0 0 0 0 0 0 0 0 0 0 0 0 0];        % 参数下限
  3. ub = [+inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf];    % 参数上限
复制代码

回复此楼 已获打赏: 0 积分

举报

新手

10 麦片

财富积分


050


1

主题

2

帖子

0

最佳答案
 楼主| 发表于 2020-7-5 13:41:30 | 显示全部楼层
20141303 发表于 2020-7-2 13:45
仅供参考,将

修改为

您好,我修改完以后会显示计算超出函数计算限值,该怎么解决呢
回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

6

主题

1440

帖子

306

最佳答案
  • 关注者: 70
发表于 2020-7-5 15:05:38 | 显示全部楼层
MATLAB2019b版本运行正常
1.PNG
回复此楼 已获打赏: 0 积分

举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /4 下一条

快速回复 返回顶部 返回列表