Fork me on GitHub

数模-梯度下降法

先占坑,太困了,明天再写


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

Adhere to original technology sharing, your support will encourage me to continue to create!