最优化计算方法

0.91
0.91
8 (x 3)2 ( y 1)2 6 (x 5)2 ( y 1)2 ] / 84
▪ 问题为在区域0=<x=<6, 0=<y=<6上求z=f(x,y)的 最小值。
绘制目标函数图形
clear all syms x y r1 = sqrt((x-1)^2+(y-5)^2)^0.91; r2 = sqrt((x-3)^2+(y-5)^2)^0.91; r3 = sqrt((x-5)^2+(y-5)^2)^0.91; r4 = sqrt((x-1)^2+(y-3)^2)^0.91; r5 = sqrt((x-3)^2+(y-3)^2)^0.91; r6 = sqrt((x-5)^2+(y-3)^2)^0.91; r7 = sqrt((x-1)^2+(y-1)^2)^0.91; r8 = sqrt((x-3)^2+(y-1)^2)^0.91; r9 = sqrt((x-5)^2+(y-1)^2)^0.91; z = 3.2+1.7*(6*r1+8*r2+8*r3+21*r4+6*r5+3*r6+18*r7+8*r8+6*r9)/84; ezmesh(z)
x1new=a+(b-a)*rand(1); x2new=c+(d-c)*rand(1); znew=subs(-z,[x1,x2],[x1new,x2new]); if znew<zmin
x1min=x1new; x2min=x2new; zmin=znew; fprintf('%4.0f %1.6f %1.6f %1.6f\n', n, x1min, x2min, zmin); end end
z f (x, y)
0.91
3.2 1.7[6 (x 1)2 ( y 5)2
0.91
0.91
8 (x 2)2 ( y 5)2 8 (x 5)2 ( y 5)2
0.91
0.91
21 (x 1)2 ( y 3)2 6 (x 3)2 ( y 3)2
0.91
0.91
3 (x 5)2 ( y 3)2 18 (x 1)2 ( y 1)2
16/5+...+17/140 (x2-10 x+26+y2-2 y)91/200
20
15
10
5
5 0
5 0
-5
-5
y
x
绘制等值线图
ezcontourf(z,[0 6 0 6])
colorbar, grid on
16/5+...+17/140 (x2-10 x+26+y2-2 y)91/200 6
牛顿法求较精确的近似值
▪ 牛顿法见书p.56 x=[x1;x2]; F=diff(z,x1); G=diff(z,x2); Dz=[F;G]; % the gradient vector of z ▪ 即求方程组Dz=0的解。
牛顿法代码实现
dFdx1=diff(F,x1); dFdx2=diff(F,x2); dGdx1=diff(G,x1); dGdx2=diff(G,x2); D2z =[dFdx1 dFdx2; dGdx1 dGdx2]; % Jacobian of Dz (same as Hessian of
zmin=f(x,y) 对n=1到N循环
开始
x=random{[a,b]}
y=random{[c,d]}
z=f(x,y) 若z<zmin,则 xmin=x,ymin=y,zmin=z 结束 结束 输出：xmin,ymin,zmin
代码实现
a=0; b=6; c=0; d=6; N=1000; x0 = a+(b-a)*rand(1); y0 = c+(d-c)*rand(1); zmin = subs(z,[x,y],[x0,y0]); fprintf(' Iteration xmin ymin zmin value\n\n'); for n=1:N
135
130
125
120
0
5
10
15
20
25
30
35
40
x
ezplot(y,[18,22]); grid on
139.4
(130-2 x) exp(1/40 x)-9/20 x
139.35
139.3
139.25
139.2
139.15
18 18.5 19 19.5 20 20.5 21 21.5 22 x
ezplot(y,[0,20])
(130-2 x) exp(1/40 x)-9/20 x 140 139 138 137 136 135 134 133 132 131 130
0
2
4
6
8 10 12 14 16 18 20
x
ezplot(y,[0,40])
(130-2 x) exp(1/40 x)-9/20 x 140
据的统计分析给出：对离救火站r英里打来
的求救电话，需要的响应时间估计
为
。下图给出了从消3.防2 1管.7r0员.91 处得到
的从城区不同区域打来的求救电话频率的
估计数据。求新的消防站的最佳位置。
301421 211232 533012 852100 10 6 3 1 3 1 023111
▪ 设(x,y)为新消防站的位置，对求救电话的平均响 应时间为：
ezplot(y,[19,20]); grid on
(130-2 x) exp(1/40 x)-9/20 x 139.395 139.394 139.393 139.392 139.391 139.39 139.389 139.388 139.387 139.386 139.385
19 19.1 19.2 19.3 19.4 19.5 19.6 19.7 19.8 19.9 20 x
xnew=a+(b-a)*rand(1); ynew=c+(d-c)*rand(1); znew=subs(z,[x,y],[xnew,ynew]); if znew<zmin
xmin=xnew; ymin=ynew; zmin=znew; fprintf('%4.0f %1.6f %1.6f %1.6f\n', n, xmin, ymin, zmin); end end
xnew=a+(b-a)*rand(1); ynew=c+(d-c)*rand(1); znew=subs(z,[x,y],[xnew,ynew]); if znew<zmin
xmin=xnew; ymin=ynew; zmin=znew; fprintf('%4.0f %1.6f %1.6f %1.6f\n', n, xmin, ymin, zmin); end end
F = dydx; F1 = diff(F,x); format long N = 10; % number of iterations x0 = 19 % initial guess fprintf(' iteration xvalue\n\n'); for i=1:N x1=x0-subs(F,x,x0)/subs(F1,x,x0); fprintf('%5.0f %1.16f\n', i, x1); x0 = x1; end display('Hence, the critical point (solution of F=0) is
(approx)'), x1
灵敏性分析
▪ 考虑最优售猪时间关于小猪增长率c=0.025 的灵敏性。
xvalues = 0; for c = 0.022:0.001:0.028 y = (0.65-0.01*x)*200*exp(c*x)-0.45*x; dydx=diff(y,x); xmaxc=solve(dydx); xmaxc = double(xmaxc); xmaxc = xmaxc(1); xvalues = [xvalues; xmaxc]; end xvalues = xvalues(2:end);
数值方法求解--Matlab
dydx = diff(y,x) xmax = solve(dydx); xmax = double(xmax) xmax =xmax(1) ymax=subs(y,x,xmax)
Newton 法
▪ 求方程F(x)=0的根. ▪ 牛顿法： x(n)=x(n-1)-F(x(n-1))/F’(x(n-1))
10 10
8 6 4 2 x2
10 8 6 4 2
x1
随机搜索求近似最优值
a=0; b=10; c=0; d=10; N=1000; x10 = a+(b-a)*rand(1); x20 = c+(d-c)*rand(1); zmin = subs(-z,[x1,x2],[x10,x20]); fprintf(' Iteration x1min x2min zmin value\n\n'); for n=1:N
例3.3 一家草坪家俱厂商生产两种草坪椅。 一种是木架的，一种是铝管架的。木架椅 的生产价格为每把18美元，铝管椅为每把 10美元。在产品出售的市场上，可以售出 的数量依赖于价格。据估计，若每天售出x 把木架椅，y把铝管椅，木架椅和铝管椅的 出售价格分别不能超过
10 31x0.5 1.3y0.2 ,5 15y0.4 0.8x0.0，8 求最优生产 量。
18*x1+ x2*(5+15*x2^(-0.4)+.8*x1^(-0.08))10*x2; ezsurfc(z, [0.1 10 0.1 10]); title('Objective Function z');