应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 求Matlab指导
101/1返回列表
查看: 809  |  回复: 9
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

zhaoshazhu

新虫 (小有名气)

[求助] 求Matlab指导 已有2人参与

为什么有时候Matlab给初值后有时计算快有时计算很慢，一直busy状态？是出现死循环还是程序有问题？请求帮忙。
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2014-12-06 21:43:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

戴钢盔的猪头

木虫 (知名作家)

【答案】应助回帖

感谢参与，应助指数 +1
初值给的不好，迭代点的二阶导数矩阵近似奇异，求解费时。

[ 发自手机版 http://muchong.com/3g ]
赞一下 回复此楼
2楼2014-12-06 22:01:04
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

yongcailiu

金虫 (小有名气)

【答案】应助回帖

感谢参与，应助指数 +1
迭代法的收敛速度及结果依赖于初值的选取，比如牛顿法，它严重依赖于初值的选取，所以初值的恰当选取对于迭代法来说还是挺重要的
赞一下 回复此楼
3楼2014-12-07 09:09:11
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhaoshazhu

新虫 (小有名气)
引用回帖:
3楼: Originally posted by yongcailiu at 2014-12-07 09:09:11
迭代法的收敛速度及结果依赖于初值的选取，比如牛顿法，它严重依赖于初值的选取，所以初值的恰当选取对于迭代法来说还是挺重要的

谢谢您的指导，我现在就是初值选取有困难，很多时候都是自己尝试，但是计算的结果不好，所得到参数的置信区间很大，请问您有好的建议吗？下面是我做动力学参数所编程的代码。可以帮我看一下有没有问题吗?
function DMMwan180
clear all
clc
k0 = [1 1 1 1 3 1 1 1 0.3 0.007];
lb = [0 0 0 0 0 0 0 0 0 0];
ub = [inf  inf  inf  inf inf inf inf inf inf inf];

P0 =[0.015481         4.644213         0.340306         0        0        0;
     0.015481         4.644213         0.340306         0        0        0;
     0.018316         4.579047         0.402637         0        0        0;
     0.022423         4.484655         0.492922         0        0        0;
     0.018316         4.579047         0.402637         0        0        0;
     0.022423         4.484655         0.492922         0        0        0;
     0.040656         4.065616         0.893728         0        0        0;
     0.015481         4.644213         0.340306         0        0        0;
     0.015481         4.644213         0.340306         0        0        0;
     0.018316         4.579047         0.402637         0        0        0;
     0.022423         4.484655         0.492922         0        0        0;
     0.040656         4.065616         0.893728         0        0        0
];  % 初始分压,MPa

Pi=[0.0000439         4.6409577         0.3473380         0.0001568         0.0058211         0.0056825;
    0.0000709         4.6414433         0.3484888         0.0003922         0.0067384         0.0028664;
    0.0000841         4.5751124         0.4132686         0.0005222         0.0086355         0.0023772;
    0.0002438         4.4758947         0.5083807         0.0014471         0.0118155         0.0022182;
    0.0003697         4.5710387         0.4156405         0.0020687         0.0094840         0.0013984;
    0.0009699         4.4759239         0.5083459         0.0032029         0.0100967         0.0014606;
    0.0036297         4.0578876         0.9185219         0.0065201         0.0116013         0.0018393;
    0.0001993         4.6368235         0.3512604         0.0010305         0.0088804         0.0018058;
    0.0009229         4.6361031         0.3516653         0.0023232         0.0080420         0.0009435;
    0.0015467         4.5720515         0.4148248         0.0039501         0.0068821         0.0007448;
    0.0036621         4.4734942         0.5082864         0.0048581         0.0086486         0.0010507;
    0.0087555         4.0513963         0.9187905         0.0087188         0.0105949         0.0017440
];
% 经过Wc/F0后，各物质分压,MPa

% 使用函数lsqnonlin()进行参数估计
opt=optimset('Algorithm','levenberg-marquardt');
[k,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(@ObjFunc,k0,[],[],opt,P0,Pi);      
ci = nlparci(k,residual,jacobian);

fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
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('\t残差平方和 = %.4f\n',resnorm)
fprintf('\texitflag = %.4f\n',exitflag)
fprintf('\tresidual = %.4f\n',residual)


% ------------------------------------------------------------------
function f = ObjFunc(k,P0,Pi)           % 目标函数
[m,n] = size(P0);
Pcal = zeros(m,n);
tspan = [0  2.05;
         0  1.37;
         0  1.62;
         0  1.98;
         0  1.21;
         0  1.49;
         0  2.70;
         0  1.02;
         0  0.82;
         0  0.97;
         0  1.19;
         0  2.16
];         % 即Wc/F0,g.h/mol
for i = 1:m
[t PP] = ode45(@Euqations,tspan(i,,P0(i,,[],k);
Pcal(i, = PP(end,;
end
f= Pcal-Pi;

% ------------------------------------------------------------------
function dPdt = Euqations(t, P, k)        % here t = Wc / F0
denom = 1+k(4)*P(1)+k(5)*P(3)+k(6)*P(4)+k(7)*P(5)+k(8)*P(6);               % k(4) = KDMM,k(5) = KME ,k(6)=KHPM,k(7)=KPDO,k(8)=KNPA,k(9)=Kp1,k(10)=Kp2
theA =k(4)*(P(1)*P(2)-P(4)*P(3)/k(9)*P(2)) / denom^2;
theB =k(6)*( P(4)*P(2)-P(5)*P(3)/k(10)*P(2))/ denom^2;
theC =k(7)*P(5)*P(2)/denom^2;
r1 = k(1)*theA;
r2 = k(2)*theB;
r3 = k(3)*theC;


dPDMMdt = -r1;
dPHdt = -2*r1-2*r2-r3;   
dPMEdt = r1+r2;
dPHPMdt = r1-r2;
dPPDOdt = r2-r3;
dPNPAdt = r3;

dPdt = [dPDMMdt;dPHdt;dPMEdt;dPHPMdt;dPPDOdt;dPNPAdt];
赞一下 回复此楼
4楼2014-12-07 10:36:02
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhaoshazhu

新虫 (小有名气)
引用回帖:
2楼: Originally posted by 戴钢盔的猪头 at 2014-12-06 22:01:04
初值给的不好，迭代点的二阶导数矩阵近似奇异，求解费时。

请问你对初值的选取有什么好的指导建议吗？因为我没有什么参考文献好多时候都是自己尝试，模拟出来的效果不好。
赞一下 回复此楼
5楼2014-12-07 10:38:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

戴钢盔的猪头

木虫 (知名作家)

【答案】应助回帖

引用回帖:
5楼: Originally posted by zhaoshazhu at 2014-12-07 10:38:09
请问你对初值的选取有什么好的指导建议吗？因为我没有什么参考文献好多时候都是自己尝试，模拟出来的效果不好。...

不同的问题不同对待，凸规划不必纠结初值，必定收敛，其它情况，如果目标函数能够通过一定的操作转化成线性问题，可以用相应线性问题的解作为初值，其它实在无法简化的非线性问题，可以在优化前用全局搜索算法搜索一个近似解，再以近似解作为初值。

[ 发自手机版 http://muchong.com/3g ]
赞一下 回复此楼
6楼2014-12-07 14:54:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhaoshazhu

新虫 (小有名气)
引用回帖:
6楼: Originally posted by 戴钢盔的猪头 at 2014-12-07 14:54:12
不同的问题不同对待，凸规划不必纠结初值，必定收敛，其它情况，如果目标函数能够通过一定的操作转化成线性问题，可以用相应线性问题的解作为初值，其它实在无法简化的非线性问题，可以在优化前用全局搜索算法搜索 ...

请问怎么用全局搜索算法找近似值啊？我对Matlab只知道皮毛，这是我毕业论文动力学模拟的一部分，所以希望得到您的帮助。我算了很长时间都没有结果，很着急。
赞一下 回复此楼
7楼2014-12-07 14:58:46
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhaoshazhu

新虫 (小有名气)
引用回帖:
6楼: Originally posted by 戴钢盔的猪头 at 2014-12-07 14:54:12
不同的问题不同对待，凸规划不必纠结初值，必定收敛，其它情况，如果目标函数能够通过一定的操作转化成线性问题，可以用相应线性问题的解作为初值，其它实在无法简化的非线性问题，可以在优化前用全局搜索算法搜索 ...

这是我的算法您可以参考一下，希望得到您的帮助。
function DMMwan180
clear all
clc
k0 = [10 1 1 1 3 10 1 1 0.3 0.007];
lb = [0 0 0 0 0 0 0 0 0 0];
ub = [inf  inf  inf  inf inf inf inf inf inf inf];

P0 =[0.015481         4.644213         0.340306         0        0        0;
     0.015481         4.644213         0.340306         0        0        0;
     0.018316         4.579047         0.402637         0        0        0;
     0.022423         4.484655         0.492922         0        0        0;
     0.018316         4.579047         0.402637         0        0        0;
     0.022423         4.484655         0.492922         0        0        0;
     0.040656         4.065616         0.893728         0        0        0;
     0.015481         4.644213         0.340306         0        0        0;
     0.015481         4.644213         0.340306         0        0        0;
     0.018316         4.579047         0.402637         0        0        0;
     0.022423         4.484655         0.492922         0        0        0;
     0.040656         4.065616         0.893728         0        0        0
];  % 初始分压,MPa

Pi=[0.0000439         4.6409577         0.3473380         0.0001568         0.0058211         0.0056825;
    0.0000709         4.6414433         0.3484888         0.0003922         0.0067384         0.0028664;
    0.0000841         4.5751124         0.4132686         0.0005222         0.0086355         0.0023772;
    0.0002438         4.4758947         0.5083807         0.0014471         0.0118155         0.0022182;
    0.0003697         4.5710387         0.4156405         0.0020687         0.0094840         0.0013984;
    0.0009699         4.4759239         0.5083459         0.0032029         0.0100967         0.0014606;
    0.0036297         4.0578876         0.9185219         0.0065201         0.0116013         0.0018393;
    0.0001993         4.6368235         0.3512604         0.0010305         0.0088804         0.0018058;
    0.0009229         4.6361031         0.3516653         0.0023232         0.0080420         0.0009435;
    0.0015467         4.5720515         0.4148248         0.0039501         0.0068821         0.0007448;
    0.0036621         4.4734942         0.5082864         0.0048581         0.0086486         0.0010507;
    0.0087555         4.0513963         0.9187905         0.0087188         0.0105949         0.0017440
];
% 经过Wc/F0后，各物质分压,MPa

% 使用函数lsqnonlin()进行参数估计
opt=optimset('Algorithm','levenberg-marquardt');
[k,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(@ObjFunc,k0,[],[],opt,P0,Pi);      
ci = nlparci(k,residual,jacobian);

fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
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('\t残差平方和 = %.4f\n',resnorm)
fprintf('\texitflag = %.4f\n',exitflag)
fprintf('\tresidual = %.4f\n',residual)


% ------------------------------------------------------------------
function f = ObjFunc(k,P0,Pi)           % 目标函数
[m,n] = size(P0);
Pcal = zeros(m,n);
tspan = [0  2.05;
         0  1.37;
         0  1.62;
         0  1.98;
         0  1.21;
         0  1.49;
         0  2.70;
         0  1.02;
         0  0.82;
         0  0.97;
         0  1.19;
         0  2.16
];         % 即Wc/F0,g.h/mol
for i = 1:m
[t PP] = ode45(@Euqations,tspan(i,,P0(i,,[],k);
Pcal(i, = PP(end,;
end
f= Pcal-Pi;

% ------------------------------------------------------------------
function dPdt = Euqations(t, P, k)        % here t = Wc / F0
denom = 1+k(4)*P(1)+k(5)*P(3)+k(6)*P(4)+k(7)*P(5)+k(8)*P(6);               % k(4) = KDMM,k(5) = KME ,k(6)=KHPM,k(7)=KPDO,k(8)=KNPA,k(9)=Kp1,k(10)=Kp2
theA =k(4)*(P(1)*P(2)-P(4)*P(3)/k(9)*P(2)) / denom^2;
theB =k(6)*( P(4)*P(2)-P(5)*P(3)/k(10)*P(2))/ denom^2;
theC =k(7)*P(5)*P(2)/denom^2;
r1 = k(1)*theA;
r2 = k(2)*theB;
r3 = k(3)*theC;


dPDMMdt = -r1;
dPHdt = -2*r1-2*r2-r3;   
dPMEdt = r1+r2;
dPHPMdt = r1-r2;
dPPDOdt = r2-r3;
dPNPAdt = r3;

dPdt = [dPDMMdt;dPHdt;dPMEdt;dPHPMdt;dPPDOdt;dPNPAdt];
赞一下 回复此楼
8楼2014-12-07 15:03:22
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

戴钢盔的猪头

木虫 (知名作家)

【答案】应助回帖

引用回帖:
7楼: Originally posted by zhaoshazhu at 2014-12-07 14:58:46
请问怎么用全局搜索算法找近似值啊？我对Matlab只知道皮毛，这是我毕业论文动力学模拟的一部分，所以希望得到您的帮助。我算了很长时间都没有结果，很着急。...

MATLAB吗？你随便用一种全局搜索算法，比如遗传算法，模拟退火或者粒子群等，具体命令请自行百度，如果所有未知量都在比较接近的数量级就随便设个初值开始全局搜索，然后把结果作为后面经典优化的初值。如果未知量之间数量级差别较大，可将未知量x改写为10^y，便于在数量级跨度大的情况下进行全局搜索。

[ 发自手机版 http://muchong.com/3g ]
赞一下 回复此楼
9楼2014-12-07 15:07:23
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhaoshazhu

新虫 (小有名气)
引用回帖:
9楼: Originally posted by 戴钢盔的猪头 at 2014-12-07 15:07:23
MATLAB吗？你随便用一种全局搜索算法，比如遗传算法，模拟退火或者粒子群等，具体命令请自行百度，如果所有未知量都在比较接近的数量级就随便设个初值开始全局搜索，然后把结果作为后面经典优化的初值。如果未知量 ...

谢谢您的指导，遗传算法表示看不懂。
赞一下 回复此楼
10楼2014-12-07 16:10:43
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 zhaoshazhu 的主题更新
101/1返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭