蒙特卡罗方法求解有约束的非线性规划问题的matlab程序
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
蒙特卡罗方法求解有约束的非线性规划问题的matlab程序
首先我们说明一下：观察函数f(x)，它恰好可以应用雅克比矩阵来计算得到：
clear;
syms x1 x2 x3 x4 x5;
y=(1+(x1*x2+x3+x4^2)^(1/2)/(x5+x3^2))^(1/2);
f=simple(jacobian(y)*[x1;x2;x3;x4;x5]);
然后将得到的f表达式前面加上@(x1,x2,x3,x4,x5)：
f=@(x1,x2,x3,x4,x5)-1/4*x3*(2*x1*x2*x3+x5+3*x3^2+2*x4^2*x3)/(x5+x3^2)/((x5+x3^2+(x1 *x2+x3+x4^2)^(1/2))*(x5+x3^2))^(1/2)/(x1*x2+x3+x4^2)^(1/2);
作为下面程序中的f即可，已经验证过结果仍然差不多。
下面程序中的f仍为原程序所给的f。
matlab蒙特卡罗方法程序（相关原理参见附件文献）：
clear;
f=@(x1,x2,x3,x4,x5)-(x3*(3*x3^2 + 2*x3*x4^2 + 2*x1*x2*x3 + x5))...
/(4*((x5 + (x4^2 + x3 + x1*x2)^(1/2) + x3^2)/(x3^2 + x5))^(1/2)*(x3^2 + x5)^2*(x4^2 + x3 + x1*x2)^(1/2));
MIN=inf;
LIMIT=10000;
while LIMIT>0
x(3)=4.5*rand+0.5;%将【0,1】区间上的随机数转化到【0.5,5】上的随机数，下面其余数类同
x(4)=2*rand+1;
x(5)=3*rand+1;
x(1)=4.5*rand+0.5;
x(2)=4.5*rand+0.5;
while x(1)+x(2)^2>=1&x(1)+x(2)^2<=10
x1=x(1);x2=x(2);x3=x(3);x4=x(4);x5=x(5);
temp=f(x1,x2,x3,x4,x5);
if temp MIN=temp; x0=x; end; break; end; LIMIT=LIMIT-1; end; MIN,x0 或者 clear; f=@(x1,x2,x3,x4,x5)-(x3*(3*x3^2 + 2*x3*x4^2 + 2*x1*x2*x3 + x5))... /(4*((x5 + (x4^2 + x3 + x1*x2)^(1/2) + x3^2)/(x3^2 + x5))^(1/2)*(x3^2 + x5)^2*(x4^2 + x3 + x1*x2)^(1/2)); MIN=inf; LIMIT=10000; while LIMIT>0 x(3)=4.5*rand+0.5; x(4)=2*rand+1; x(5)=3*rand+1; x(1)=4.5*rand+0.5; x(2)=4.5*rand+0.5; if x(1)+x(2)^2>=1&x(1)+x(2)^2<=10 x1=x(1);x2=x(2);x3=x(3);x4=x(4);x5=x(5); temp=f(x1,x2,x3,x4,x5); if temp MIN=temp; x0=x; end; end LIMIT=LIMIT-1; end; MIN,x0 前者几次模拟计算结果如下： MIN = -0.3003 x0 = 4.6885 2.2581 1.4226 2.4956 1.1308 MIN = -0.3135 x0 = 2.9750 2.2140 1.3397 2.8936 1.0118 MIN = -0.3173 x0 = 4.9725 2.2098 1.5285 2.7026 1.0338 注：源程序如下 clear; f=@(x1,x2,x3,x4,x5)-(x3*(3*x3^2 + 2*x3*x4^2 + 2*x1*x2*x3 + x5))... /(4*((x5 + (x4^2 + x3 + x1*x2)^(1/2) + x3^2)/(x3^2 + x5))^(1/2)*(x3^2 + x5)^2*(x4^2 + x3 + x1*x2)^(1/2)); MIN=inf; LIMIT=10000; while LIMIT>0 x(3)=4.5*rand+0.5; x(4)=2*rand+1; x(5)=3*rand+1; x(1)=4.5*rand+0.5; x(2)=4.5*rand+0.5; while x(1)+x(2)^2<1 | x(1)+x(2)^2>10 x(1)=4.5*rand+0.5; x(2)=4.5*rand+0.5; end; x1=x(1);x2=x(2);x3=x(3);x4=x(4);x5=x(5); temp=f(x1,x2,x3,x4,x5); if temp MIN=temp; x0=x; end; LIMIT=LIMIT-1;