先占坑,太困了,明天再写
2018/4/7 20:04:14 更新
感觉蛮有用的,放在这以供查看
数学建模第四次课程作业
非线性规划
用梯度下降法求解 min f (x) = x1 - x2 + 2x12 + 2x1x2 + x22,初始点取x0 = (0, 0)T,允许误差取e = 0.05。
注:只需迭代三步,即求到x(3)。
解:
手工计算到x(3)
(1) 计算目标函数f (x)的负梯度,有:
-▽f (x) = ((-1-4x1-2x2)¦(1-2x1-2x2))
(2)选择初始点x(0) = (0,0)T,计算得:
P(0) = -▽f (x(0)) = (-1,1)T,‖P0‖ = √2 ≈ 1.4 > e = 0.05
故x(0)不能成为近似最优解,需考虑最优步长λ0。
将x(0) + λP(0) = ((-λ)¦λ)代入 f (x)有f (x(0) +λP(0)) =λ2 - 2λ。求单变量极值问题得:
λ0 = (min)┬(λ≥0)〖f (x(0) + λP(0) 〗) =1,于是由迭代格式有x(1) = x(0) + λP(0) = ((-1)¦1)
(3)对于点x(1) = (-1,1)T,计算得:
P(1) = -▽f (x(1)) = (1,1)T,‖P1‖ = √2 ≈ 1.4 > e = 0.05
故x(1)不能成为近似最优解,需考虑最优步长λ1。
将x(1) +λP(1) = ((-1+λ)¦(1+λ))代入 f (x)有f (x(1) +λP(1)) = 5λ2 - 2λ - 1。求单变量极值问题得:
λ1 = (min)┬(λ≥0)〖f (x(1) +λP(1)) 〗 =1/5,于是由迭代格式有x(2) = x(1) +λP(1) = ((-4/5)¦(6/5))
(4)对于点x(2) = (-4/5,6/5)T,计算得:
P(2) = -▽f (x(2)) = (-1/5,1/5)T,‖P2‖ = √2/5 ≈ 0.28 > e = 0.05
故x(2)不能成为近似最优解,需考虑最优步长λ2。
将x(2) +λP(2) = ((-4/5-1/5 λ)¦(6/5+1/5 λ))代入 f (x)有f (x(2) +λP(2)) = 1/25 (λ-1)2 - 31/25 。求单变量极值问题得:
λ2 = (min)┬(λ≥0)〖f (x(2) +λP(2)) 〗 = 1,于是由迭代格式有x(3) = x(2) +λP(2) = ((-1)¦(7/5))
(5)对于点x(3) = (-1,7/5)T,计算得:
P(3) = -▽f (x(3)) = (1/5,1/5)T,‖P3‖ = √2/5 ≈ 0.28 > e = 0.05
故x(3)不能成为近似最优解,需考虑最优步长λ2。
将x(2) +λP(2) = ((-1+1/5 λ)¦(7/5+1/5 λ))代入 f (x)有f (x(2) +λP(2)) = 1/5 (λ- 6/5)2 - 31/25 。求单变量极值问题得:
λ3 = (min)┬(λ≥0)〖f (x(3) +λP(3)) 〗 = 6/5,于是由迭代格式有x(4) = x(3) +λP(3) = ((-19/25)¦(41/25))
(6) ······
Matlab函数实现【参考博客】:
函数部分
function [ y ] = steepest(fx,var,x,e,MAX)
if nargin < 5
MAX = 10;
end
precision = 3;
bfound = 0;
for k=1:1:MAX
direction = getNextDirecrion(fx,var,x);
disp('------------------------------');
fprintf('d[%d]=:',k);
disp( vpa(direction',precision) );
if normest(direction) <= e
y = x;
bfound = 1;
break;
else
step = getNextStep(fx,var, x,direction);
if isempty(step)
error('can not find a proper step.');
end
%´òÓ¡Çó½â¹ý³Ì
fprintf('X[%d]=:',k);
disp( vpa(x',precision) );
fprintf('step(%d)=: ', k);
disp( vpa(step,precision) );
disp('------------------------------');
x = x+step*direction;
end
end
if bfound == 1
disp('min value of:');
disp( vpa( subs(fx,var,y),precision) );
end
end
function [direction] = getNextDirecrion(fx,var,xk)
gx = gradient(fx,var);
direction = -subs(gx,var,xk);
end
function [step] =getNextStep(fx,var,xk,dk)
syms lambda;
phix = subs(fx,var,xk+lambda*dk);
phix_diff = diff(phix);
step = double(solve(phix_diff,'Real',true));
end
Input:
>> syms x1 x2;
X = [x1;x2];
fx = x1-x2+2*x1^2+2*x1*x2+x2^2;
x1 = [0;0];
e = 0.05;
>> minVal = steepest(fx,X,x1,e)
Output:
------------------------------
d[1]=:[ -1.0, 1.0]
X[1]=:[ 0, 0]
step(1)=: 1.0
------------------------------
------------------------------
d[2]=:[ 1.0, 1.0]
X[2]=:[ -1.0, 1.0]
step(2)=: 0.2
------------------------------
------------------------------
d[3]=:[ -0.2, 0.2]
X[3]=:[ -0.8, 1.2]
step(3)=: 1.0
------------------------------
------------------------------
d[4]=:[ 0.2, 0.2]
X[4]=:[ -1.0, 1.4]
step(4)=: 0.2
------------------------------
------------------------------
d[5]=:[ -0.04, 0.04]
X[5]=:[ -0.96, 1.44]
step(5)=: 1.0
------------------------------
------------------------------
d[6]=:[ 0.04, 0.04]
X[6]=:[ -1.0, 1.48]
step(6)=: 0.2
------------------------------
------------------------------
d[7]=:[ -0.008, 0.008]
min value of:
-1.25
minVal =
-124/125
186/125
结论
故需要迭代7步得到最优值及最优解,最优值为x = (-124/125,186/125),最优解为f (x) = -1.25