应用最优化方法及MATLAB实现——第4章代码实现

一、概述

        之前对这本书的第三章进行了代码实现，这篇博客是对这本书第4章相关代码进行实现，部分内容安装书中代码无法实现相应功能，MATLAB会报错，对其进行一定程度的更改后，可以正常运行，与书中所给示例运行结果相一致。

二、具体实现

（一）使用函数

        1.概述

        因为书中几个条件使用的示例函数均相同，在此将其单独罗列出来，复制即可。

        2.函数实现

        f_test1和其导数g_test1。

function y = f_test1(x)
%F_TEST1 此处显示有关此函数的摘要
%   此处显示详细说明y = -3 * x * sin(0.75 * x) + exp(-2 * x);end

function y = g_test1(x)
%G_TEST1 此处显示有关此函数的摘要
%   此处显示详细说明y = -2 / exp(2 * x) - 3 * sin(x * 0.75) - (3 * 0.75 * x * cos(0.75 * x));end

        f_test2和其导数g_test2。

function y = f_test2(x)
%F_TEST2 此处显示有关此函数的摘要
%   此处显示详细说明x1 = x(1);
x2 = x(2);
y = x1 ^ 2 + x2 ^ 2 - 1;end

function y = g_test2(x)
%G_TEST2 此处显示有关此函数的摘要
%   此处显示详细说明x1 = x(1);
x2 = x(2);
y = 2 * x1 + 2 * x2;end

         f_test3和其导数g_test3。

function y = f_test3(x)
%F_TEST3 此处显示有关此函数的摘要
%   此处显示详细说明y = sin(3 * x) / x;end

function y = g_test3(x)
%G_TEST3 此处显示有关此函数的摘要
%   此处显示详细说明y = (3 * cos(3 * x)) / x - sin(3 * x) / ( x ^ 2);end

（二）Armijo_search条件

        1.main.m文件

        这个文件是此条件的主运行文件，放开相应注释即可运行每个示例。

% 这个文件主要为Armijo_search文件的主程序% 清空
close;
clear;
clc;% 第一个示例
% x_current = -2;
% d_current = 1;
% rho = 0.1;
% [alpha_acceptable, x_next, f_next, k] = Armijo_search(@f_test1, @g_test1, x_current, d_current, rho);% 第二个示例
% x_current = [2;2];
% d_current = [-1;-1];
% rho = 0.1;
% [alpha_acceptable, x_next, f_next, k] = Armijo_search(@f_test2, @g_test2, x_current, d_current, rho);% 第三个示例
x_current = 5;
d_current = 1;
rho = 0.5;
[alpha_acceptable, x_next, f_next, k] = Armijo_search(@f_test3, @g_test3, x_current, d_current, rho);

        2.Armijo_search.m文件

        此函数跟书中函数一样，没有太大改动。

function [alpha_acceptable, x_next, f_next, k] = Armijo_search(f_test, g_test, x_current, d_current, rho)
% f_test, 目标函数
% g_test, 目标函数对决策变量x的导数
% x_current, x在向量空间中的当前点
% d_current, f_test在x_current的下降搜索方向
% rho, 可接受系数k_max = 1000;
k = 0;
f_current = f_test(x_current);
g_current = g_test(x_current);
f_alpha_lower_k = f_current;
g_alpha_lower_k = g_current;
df_alpha_lower_k = (d_current') * g_alpha_lower_k; % 这里这个值是一直保持不变的
f_alpha_lower_0 = f_alpha_lower_k;
df_alpha_lower_0 = df_alpha_lower_k;alpha_lower_k = 0;
alpha_upper_k = 1e8;
alpha_k = alpha_upper_k;% 
for k = 1:k_maxx_alpha_k = x_current + alpha_k * d_current;f_alpha_k = f_test(x_alpha_k);Armijo_condition = f_alpha_k - f_alpha_lower_0 - rho * alpha_k * df_alpha_lower_0;if (Armijo_condition <= 0)alpha_acceptable = alpha_k;x_next = x_alpha_k;f_next = f_alpha_k;break;elseif (alpha_k < alpha_upper_k)alpha_upper_k = alpha_k;endalpha_k = alpha_lower_k + (1/2) * ((alpha_k ^ 2) * df_alpha_lower_k) / (f_alpha_lower_k - f_alpha_k + alpha_k * df_alpha_lower_k);% x_alpha_k = x_current + alpha_k * d_current;% g_alpha_lower_k = g_test(x_alpha_k);% df_alpha_lower_k = (d_current') * g_alpha_lower_k;end
endif(k == k_max)disp('Armijo inexact line search algorithm failed');alpha_acceptable = NaN;x_next = NaN;f_next = NaN;
endend

（三）Goldstein_search条件

        1.main.m文件

        这个文件是此条件的主运行文件，放开相应注释即可运行每个示例。

% 这个文件主要为Goldstein_search文件的主程序% 清空
close;
clear;
clc;% 第一个示例
% x_current = -2;
% d_current = 1;
% rho = 0.1;
% [alpha_acceptable, x_next, f_next, k] = Goldstein_search(@f_test1, @g_test1, x_current, d_current, rho);% 第二个示例
x_current = [2;2];
d_current = [-1;-1];
rho = 0.1;
[alpha_acceptable, x_next, f_next, k] = Goldstein_search(@f_test2, @g_test2, x_current, d_current, rho);% 第三个示例
% x_current = 5;
% d_current = 1;
% rho = 0.1;
% [alpha_acceptable, x_next, f_next, k] = Goldstein_search(@f_test3, @g_test3, x_current, d_current, rho);

        2.Goldstein_search.m文件

        此文件有些改动。

function [alpha_acceptable, x_next, f_next, k] = Goldstein_search(f_test, g_test, x_current, d_current, rho)
% f_test, 目标函数
% g_test, 目标函数对决策变量x的导数
% x_current, x在向量空间中的当前点
% d_current, f_test在x_current的下降搜索方向
% rho, 可接受系数k_max = 1000;
k = 0;
f_current = f_test(x_current);
g_current = g_test(x_current);
f_alpha_lower_k = f_current;
g_alpha_lower_k = g_current;
df_alpha_lower_k = (d_current') * g_alpha_lower_k; 
f_alpha_lower_0 = f_alpha_lower_k;
df_alpha_lower_0 = df_alpha_lower_k;tolerance = 1e-15;
if (abs(df_alpha_lower_k) > tolerance)alpha_initial = - 2 * f_alpha_lower_k ./ df_alpha_lower_k;
elsealpha_initial = 1;
end
alpha_lower_k = 0;
alpha_upper_k = 1e8;
alpha_k = alpha_initial; % 这个值是从初始值开始for k = 1:k_maxx_alpha_k = x_current + alpha_k .* d_current;f_alpha_k = f_test(x_alpha_k);g_alpha_k = g_test(x_alpha_k);df_alpha_k = (d_current') * g_alpha_k;Goldstein_condition1 = f_alpha_k - f_alpha_lower_0 - rho * alpha_k * (df_alpha_lower_0');Goldstein_condition2 = f_alpha_lower_0 + (1 - rho) * alpha_k * (df_alpha_lower_0') - f_alpha_k;if(Goldstein_condition1 <= 0)if(Goldstein_condition2 <= 0)alpha_acceptable = alpha_k;x_next = x_alpha_k;f_next = f_alpha_k;break;elsedelta_alpha_k = (alpha_k - alpha_lower_k) * df_alpha_k / (df_alpha_lower_k - df_alpha_k);if(delta_alpha_k <= 0)alpha_k_temp = 2 * alpha_k;elsealpha_k_temp = alpha_k + delta_alpha_k;endalpha_lower_k = alpha_k;f_alpha_lower_k = f_alpha_k;df_alpha_lower_k = df_alpha_k;alpha_k = alpha_k_temp;endelseif (alpha_k < alpha_upper_k)alpha_upper_k = alpha_k;endalpha_k_temp = alpha_lower_k + (1/2) * (((alpha_k - alpha_lower_k) ^ 2) * df_alpha_lower_k) / (f_alpha_lower_k - f_alpha_k + (alpha_k - alpha_lower_k) * df_alpha_lower_k);alpha_k = alpha_k_temp;endif(alpha_upper_k - alpha_lower_k < tolerance)alpha_acceptable = alpha_k;x_next = x_alpha_k;f_next = f_alpha_k;break;end
end
if((Goldstein_condition1 > 0)||(Goldstein_condition2 > 0))disp('Goldstein inexact line search algorithm failed');alpha_acceptable = NaN;x_next = NaN;f_next = NaN;
end
end

        3.注意

        （1）第1方面        

这个文件改动的主要原因是因为，在运行第二个示例中，无法安装成功运行，会出现数据维度不一样的报错。因为第二个示例所给的数据不再是单一变量，而是双变量。

        如图所示，如图中红线划出部分。这个部分在单一变量中不会出错，但是在多变量中，因为向量没有除法，所以会报错。

        将其改成如图所示，即可消除报错。

        （2）第二方面       

        除此之外，还有一部分，如图中红线圈出，由于这两个向量维度相同，都是2*1，所以无法直接相乘。

        将其更改为如图所示的部分，更改的思路是因为最后需要的是标量，需要将后面的维度为2*1的，进行转置即可。

 （四）Wolfe_search条件

        1.main.m文件

        这个文件是此条件的主运行文件，放开相应注释即可运行每个示例。

% 这个文件主要为Wolfe_search文件的主程序% 清空
close;
clear;
clc;% 第一个示例
% x_current = -2;
% d_current = 1;
% rho = 0.1;
% sigma = 0.11;
% [alpha_acceptable, x_next, f_next, k] = Wolfe_search(@f_test1, @g_test1, x_current, d_current, rho, sigma);% 第二个示例
x_current = [2;2];
d_current = [-1;-1];
rho = 0.1;
sigma = 0.11;
[alpha_acceptable, x_next, f_next, k] = Wolfe_search(@f_test2, @g_test2, x_current, d_current, rho, sigma);% 第三个示例
% x_current = 5;
% d_current = 1;
% rho = 0.5;
% sigma = 0.11;
% [alpha_acceptable, x_next, f_next, k] = Wolfe_search(@f_test3, @g_test3, x_current, d_current, rho, sigma);

        2.Wolfe_search.m文件

        此文件的部分地方同样也经过更改，部分更改的思路与前面有一些相似之处。

function [alpha_acceptable, x_next, f_next, k] = Wolfe_search(f_test, g_test, x_current, d_current, rho, sigma)
% f_test, 目标函数
% g_test, 目标函数对决策变量x的导数
% x_current, x在向量空间中的当前点
% d_current, f_test在x_current的下降搜索方向
% rho, 可接受系数
% sigma, 可接受点处的切线斜率大于初始点处切线斜率的倍数，0<rho<sigma<1k_max = 1000;
k = 0;
f_current = f_test(x_current);
g_current = g_test(x_current);
f_alpha_lower_k = f_current;
g_alpha_lower_k = g_current;
df_alpha_lower_k = (d_current') * g_alpha_lower_k; 
f_alpha_lower_0 = f_alpha_lower_k;
df_alpha_lower_0 = df_alpha_lower_k;tolerance = 1e-15;
if (abs(df_alpha_lower_k) > tolerance)alpha_initial = - 2 * f_alpha_lower_k ./ df_alpha_lower_k;
elsealpha_initial = 1;
end
alpha_lower_k = 0;
alpha_upper_k = 1e8;
alpha_k = alpha_initial; % 这个值是从初始值开始for k = 1:k_maxx_alpha_k = x_current + alpha_k .* d_current;f_alpha_k = f_test(x_alpha_k);g_alpha_k = g_test(x_alpha_k);df_alpha_k = (d_current') * g_alpha_k;Wolfe_condition1 = f_alpha_k - f_alpha_lower_0 - rho * alpha_k * (df_alpha_lower_0');Wolfe_condition2 = sum(sigma * df_alpha_lower_0 - df_alpha_k);if(Wolfe_condition1 <= 0)if(Wolfe_condition2 <= 0)alpha_acceptable = alpha_k;x_next = x_alpha_k;f_next = f_alpha_k;break;elsedelta_alpha_k = (alpha_k - alpha_lower_k) .* df_alpha_k ./ (df_alpha_lower_k - df_alpha_k);if(delta_alpha_k <= 0)alpha_k_temp = 2 * alpha_k;elsealpha_k_temp = alpha_k + delta_alpha_k;endalpha_lower_k = alpha_k;f_alpha_lower_k = f_alpha_k;df_alpha_lower_k = df_alpha_k;alpha_k = alpha_k_temp;endelseif (alpha_k < alpha_upper_k)alpha_upper_k = alpha_k;endalpha_k_temp = alpha_lower_k + (1/2) * (((alpha_k - alpha_lower_k) ^ 2) * df_alpha_lower_k) / (f_alpha_lower_k - f_alpha_k + (alpha_k - alpha_lower_k) * df_alpha_lower_k);alpha_k = alpha_k_temp;endif(alpha_upper_k - alpha_lower_k < tolerance)alpha_acceptable = alpha_k;x_next = x_alpha_k;f_next = f_alpha_k;break;end
end
if((Wolfe_condition1 > 0)||(Wolfe_condition2 > 0))disp('Wolfe inexact line search algorithm failed');alpha_acceptable = NaN;x_next = NaN;f_next = NaN;
end
end

        3.注意

        （1）第1方面

        与Goldstein条件相同的原因，需要将其修改为下面这样，如图所示。

        （2）第2方面

        与Goldstein条件相同的原因，需要将其修改为下面这样，如图所示。原因还是因为维度问题。

         （3）第3方面

        如书中源代码所示，红线划出部分。这个Wolfe_condition2计算结果是2*1维的向量，不是一个标量，虽然在下面判断Wolfe_condition2是否≤0时候没有影响。

        但其会影响到这个部分，最后的部分的判断会报错，如图所示。

         根据对书中内容的分析，此处应该是个标量，而不应该是个向量，如图所示。

         使用sum将其每项相加即可。

（五）Strong_Wolfe_search条件

         1.main.m文件

        这个文件是此条件的主运行文件，放开相应注释即可运行每个示例。

% 这个文件主要为Strong_Wolfe_search文件的主程序% 清空
close;
clear;
clc;% 第一个示例
x_current = -2;
d_current = 1;
rho = 0.1;
sigma = 0.11;
[alpha_acceptable, x_next, f_next, k] = Strong_Wolfe_search(@f_test1, @g_test1, x_current, d_current, rho, sigma);% 第二个示例
% x_current = [2;2];
% d_current = [-1;-1];
% rho = 0.1;
% sigma = 0.11;
% [alpha_acceptable, x_next, f_next, k] = Strong_Wolfe_search(@f_test2, @g_test2, x_current, d_current, rho, sigma);% 第三个示例
% x_current = 5;
% d_current = 1;
% rho = 0.1;
% sigma = 0.11;
% [alpha_acceptable, x_next, f_next, k] = Strong_Wolfe_search(@f_test3, @g_test3, x_current, d_current, rho, sigma);

        2.Strong_Wolfe_search.m文件

        此文件经过部分更改。

function [alpha_acceptable, x_next, f_next, k] = Strong_Wolfe_search(f_test, g_test, x_current, d_current, rho, sigma)
% f_test, 目标函数
% g_test, 目标函数对决策变量x的导数
% x_current, x在向量空间中的当前点
% d_current, f_test在x_current的下降搜索方向
% rho, 可接受系数
% sigma, 可接受点处的切线斜率大于初始点处切线斜率的倍数，0<rho<sigma<1k_max = 1000;
k = 0;
f_current = f_test(x_current);
g_current = g_test(x_current);
f_alpha_lower_k = f_current;
g_alpha_lower_k = g_current;
df_alpha_lower_k = (d_current') * g_alpha_lower_k; 
f_alpha_lower_0 = f_alpha_lower_k;
df_alpha_lower_0 = df_alpha_lower_k;tolerance = 1e-15;
if (abs(df_alpha_lower_k) > tolerance)alpha_initial = - 2 * f_alpha_lower_k ./ df_alpha_lower_k;
elsealpha_initial = 1;
end
alpha_lower_k = 0;
alpha_upper_k = 1e8;
alpha_k = alpha_initial; % 这个值是从初始值开始for k = 1:k_maxx_alpha_k = x_current + alpha_k .* d_current;f_alpha_k = f_test(x_alpha_k);g_alpha_k = g_test(x_alpha_k);df_alpha_k = (d_current') * g_alpha_k;Strong_Wolfe_condition1 = f_alpha_k - f_alpha_lower_0 - rho * alpha_k * (df_alpha_lower_0');Strong_Wolfe_condition2 = sum(sigma * df_alpha_lower_0 + abs(df_alpha_k));% Strong_Wolfe_condition2 = sum(abs(df_alpha_k) - abs(sigma * df_alpha_lower_0));if(Strong_Wolfe_condition1 <= 0)if(Strong_Wolfe_condition2 <= 0)alpha_acceptable = alpha_k;x_next = x_alpha_k;f_next = f_alpha_k;break;elsedelta_alpha_k = (alpha_k - alpha_lower_k) .* df_alpha_k ./ (df_alpha_lower_k - df_alpha_k);if(delta_alpha_k <= 0)alpha_k_temp = 2 * alpha_k;elsealpha_k_temp = alpha_k + delta_alpha_k;endalpha_lower_k = alpha_k;f_alpha_lower_k = f_alpha_k;df_alpha_lower_k = df_alpha_k;alpha_k = alpha_k_temp;endelseif (alpha_k < alpha_upper_k)alpha_upper_k = alpha_k;endalpha_k_temp = alpha_lower_k + (1/2) * (((alpha_k - alpha_lower_k) ^ 2) * df_alpha_lower_k) / (f_alpha_lower_k - f_alpha_k + (alpha_k - alpha_lower_k) * df_alpha_lower_k);alpha_k = alpha_k_temp;endif(alpha_upper_k - alpha_lower_k < tolerance)alpha_acceptable = alpha_k;x_next = x_alpha_k;f_next = f_alpha_k;break;end
end
if((Strong_Wolfe_condition1 > 0)||(Strong_Wolfe_condition2 > 0))disp('Wolfe inexact line search algorithm failed');alpha_acceptable = NaN;x_next = NaN;f_next = NaN;
end
end

        3.注意

        （1）第1方面

        跟前面类似，如图所示。

        （2）第2方面

        跟前面类似，如图所示，这里将两个问题合并。

        （3）第3方面

        如图所示，第二个判断条件，可以换成注释掉的那种，根据书中的条件，第二种写法更加符合书中的条件，如图所示。