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

GFCGUO

于 2024-07-12 00:58:51 发布

阅读量394

点赞数 8

文章标签：算法

本文链接：https://blog.csdn.net/GFCLJY/article/details/140361993

版权

一、概述

在阅读最优方法及MATLAB实现后，想着将书中提供的代码自己手敲一遍，来提高自己对书中内容理解程度，巩固一下。

这部分内容主要针对第3章的内容，将其所有代码实现均手敲一遍，中间部分代码自己根据其公式有些许的调整，但最后运算结果与书中所给结果一致。

修改的部分我会在下面进行单独说明并标注上我这样做的原因。

二、具体代码

（一）计算函数文件

因为书中所给出的示例代码，所使用的函数都是一样的。所以将其单独罗列出来。

其中，f_test1.m如下所示。

function y = f_test1(x)
% 该函数的主要作用是编写函数，用于计算函数输出值
% x为函数输入值
% y为函数输出值
    y = 2 * x ^2 - x - 1;
end

f_test2.m如下所示。

function y = f_test2(x)
% 该函数的主要作用是编写函数，用于计算函数输出值
% x为函数输入值
% y为函数输出值
    x1 = x(1);
    x2 = x(2);
    y = x1 ^2 + x2 ^2 - 1;
end

f_test3.m如下所示。

function y = f_test3(x)
% 该函数的主要作用是编写函数，用于计算函数输出值
% x为函数输入值
% y为函数输出值
    if(x <= 2)
        y = -x + 3;
    else
        y = x / 2;
    end
end

这些文件每个方法使用的都是一样的。

（二）对分搜索法

1.对分搜索法函数文件

这个文件主要参照书中的，没有改动。

function [alpha_star, x_next, f_next, k] = Dichotomous_search(f_test, x_current, d_current, alpha_lower, alpha_upper, tolerance)
% 该函数针对对分搜索法的实现
% 输入参数说明
% f_test, 目标函数
% x_current, x在向量空间中的当前点（已确定）
% d_current, f_test在向量空间中的搜索方向（已确定）
% alpha_lower, 从x_current出发，沿着d_current得到一个步长的单谷区间，该单谷区间的初始下届
% alpha_upper, 从x_current出发，沿着d_current得到一个步长的单谷区间，该单谷区间的初始上届
% tolerance, 最终区间的要求宽度，精度不能小于扰动量
    
% 输出参数说明
% alpha_star, 完成对分搜索后输出的步长
% x_next, x_next = x_current + alpha_star * d_current;局部最优点
% f_next, 局部最优点对应的函数值
% k,迭代次数

    disturbance_quantity = 0;
    if(tolerance >= 1e-8)
        disturbance_quantity = 1e-9;
    else
        disturbance_quantity = 0.1 * tolerance;
    end
    
    % 初始值
    k = 0;
    alpha_lower_k = alpha_lower;
    alpha_upper_k = alpha_upper;
    alpha_left_k = (1/2)*(alpha_lower_k + alpha_upper_k) - disturbance_quantity;
    alpha_right_k = (1/2)*(alpha_lower_k + alpha_upper_k) + disturbance_quantity;
    x_alpha_left_k = x_current + alpha_left_k * d_current;
    x_alpha_right_k = x_current + alpha_right_k * d_current;
    f_alpha_left_k = f_test(x_alpha_left_k);
    f_alpha_right_k = f_test(x_alpha_right_k);

    % 进入循环
    while(abs(alpha_upper_k - alpha_lower_k) > tolerance)
        if(f_alpha_right_k > f_alpha_left_k)
            alpha_upper_k = alpha_right_k;
        elseif(f_alpha_right_k < f_alpha_left_k)
            alpha_lower_k = alpha_left_k;
        else
            alpha_upper_k = alpha_right_k;
            alpha_lower_k = alpha_left_k;
        end
        alpha_left_k = (1/2)*(alpha_lower_k + alpha_upper_k) - disturbance_quantity;
        alpha_right_k = (1/2)*(alpha_lower_k + alpha_upper_k) + disturbance_quantity;
        x_alpha_left_k = x_current + alpha_left_k * d_current;
        x_alpha_right_k = x_current + alpha_right_k * d_current;
        f_alpha_left_k = f_test(x_alpha_left_k);
        f_alpha_right_k = f_test(x_alpha_right_k);
        k = k + 1;
    end
    alpha_star = (alpha_right_k + alpha_left_k) / 2;
    x_next = x_current + alpha_star * d_current;
    f_next = f_test(x_next);
end

2.对应的主运行文件

放开相应的注释即可。

% 这个文件作为第3章的主文件

% 清空所有
close;
clear;
clc;

% 对分搜索法实现
% 第一个示例
% x_current = -0.5;
% d_current = 1;
% alpha_lower = 0;
% alpha_upper = 1;
% tolerance = 1e-4;
% [alpha_star, x_next, f_next, k] = Dichotomous_search(@f_test1, x_current, d_current, alpha_lower, alpha_upper, tolerance);

% 第二个示例
% x_current = [2, 2];
% d_current = [-1, -1];
% alpha_lower = 0;
% alpha_upper = 2;
% tolerance = 1e-6;
% [alpha_star, x_next, f_next, k] = Dichotomous_search(@f_test2, x_current, d_current, alpha_lower, alpha_upper, tolerance);

% 第三个示例
x_current = 3;
d_current = -1;
alpha_lower = 0;
alpha_upper = 2;
tolerance = 1e-4;
[alpha_star, x_next, f_next, k] = Dichotomous_search(@f_test3, x_current, d_current, alpha_lower, alpha_upper, tolerance);

（三）三点等间隔搜索法

1.三点等间隔搜索法

function [alpha_star, x_next, f_next, k] = Trisection_search(f_test, x_current, d_current, alpha_lower, alpha_upper, tolerance)
% 该函数针对三点等间隔的实现
% 输入参数说明
% f_test, 目标函数
% x_current, x在向量空间中的当前点（已确定）
% d_current, f_test在向量空间中的搜索方向（已确定）
% alpha_lower, 从x_current出发，沿着d_current得到一个步长的单谷区间，该单谷区间的初始下届
% alpha_upper, 从x_current出发，沿着d_current得到一个步长的单谷区间，该单谷区间的初始上届
% tolerance, 最终区间的要求宽度，精度不能小于扰动量
    
% 输出参数说明
% alpha_star, 完成对分搜索后输出的步长
% x_next, x_next = x_current + alpha_star * d_current;局部最优点
% f_next, 局部最优点对应的函数值
% k,迭代次数
    
% 初始化
k = 0;
% 这两个是每次迭代过程中的上下界
alpha_lower_k = alpha_lower;
alpha_upper_k = alpha_upper;

alpha_left_k = alpha_lower_k + (1/4) * (alpha_upper_k - alpha_lower_k);
alpha_middle_k = alpha_lower_k + (1/2) * (alpha_upper_k - alpha_lower_k);
alpha_right_k = alpha_lower_k + (3/4) * (alpha_upper_k - alpha_lower_k);

x_alpha_left_k = x_current + alpha_left_k * d_current;
x_alpha_middle_k = x_current + alpha_middle_k * d_current;
x_alpha_right_k = x_current + alpha_right_k * d_current;

f_alpha_left_k = f_test(x_alpha_left_k);
f_alpha_middle_k = f_test(x_alpha_middle_k);
f_alpha_right_k = f_test(x_alpha_right_k);

% 这部分代码书写是一种计算方法，不太精准，部分信息没有使用到
% while abs(alpha_upper_k -alpha_lower_k) > tolerance
%     if (f_alpha_left_k < f_alpha_middle_k) && (f_alpha_left_k < f_alpha_right_k)
%         % alpha_lower_k = alpha_lower;
%         alpha_upper_k = alpha_middle_k;
%     elseif (f_alpha_middle_k < f_alpha_left_k) && (f_alpha_middle_k < f_alpha_right_k)
%         alpha_lower_k = alpha_left_k;
%         alpha_upper_k = alpha_right_k;
%     else 
%         alpha_lower_k = alpha_middle_k;
%         % alpha_upper_k = alpha_upper
%     end
% k = k + 1;
% alpha_left_k = alpha_lower_k + (1/4) * (alpha_upper_k - alpha_lower_k);
% alpha_middle_k = alpha_lower_k + (1/2) * (alpha_upper_k - alpha_lower_k);
% alpha_right_k = alpha_lower_k + (3/4) * (alpha_upper_k - alpha_lower_k);
% 
% x_alpha_left_k = x_current + alpha_left_k * d_current;
% x_alpha_middle_k = x_current + alpha_middle_k * d_current;
% x_alpha_right_k = x_current + alpha_right_k * d_current;
% 
% f_alpha_left_k = f_test(x_alpha_left_k);
% f_alpha_middle_k = f_test(x_alpha_middle_k);
% f_alpha_right_k = f_test(x_alpha_right_k);
% end

% 第二种方法，可以重复利用三点等间隔搜索法的信息
while abs(alpha_upper_k -alpha_lower_k) > tolerance
    if (f_alpha_left_k < f_alpha_middle_k) && (f_alpha_left_k < f_alpha_right_k)
        % alpha_lower_k = alpha_lower;
        alpha_upper_k = alpha_middle_k;

        alpha_left_k = alpha_lower_k + (1/4) * (alpha_upper_k - alpha_lower_k);
        alpha_middle_k = alpha_left_k;
        alpha_right_k = alpha_lower_k + (3/4) * (alpha_upper_k - alpha_lower_k);

        x_alpha_left_k = x_current + alpha_left_k * d_current;
        x_alpha_middle_k = x_current + alpha_middle_k * d_current;
        x_alpha_right_k = x_current + alpha_right_k * d_current;

        f_alpha_left_k = f_test(x_alpha_left_k);
        f_alpha_middle_k = f_alpha_left_k;
        f_alpha_right_k = f_test(x_alpha_right_k);
    elseif (f_alpha_middle_k < f_alpha_left_k) && (f_alpha_middle_k < f_alpha_right_k)
        alpha_lower_k = alpha_left_k;
        alpha_upper_k = alpha_right_k;

        alpha_left_k = alpha_lower_k + (1/4) * (alpha_upper_k - alpha_lower_k);
        % alpha_middle_k = alpha_left_k;
        alpha_right_k = alpha_lower_k + (3/4) * (alpha_upper_k - alpha_lower_k);

        x_alpha_left_k = x_current + alpha_left_k * d_current;
        % x_alpha_middle_k = x_current + alpha_middle_k * d_current;
        x_alpha_right_k = x_current + alpha_right_k * d_current;

        f_alpha_left_k = f_test(x_alpha_left_k);
        % f_alpha_middle_k = f_alpha_left_k;
        f_alpha_right_k = f_test(x_alpha_right_k);
    else 
        alpha_lower_k = alpha_middle_k;
        % alpha_upper_k = alpha_upper

        alpha_left_k = alpha_lower_k + (1/4) * (alpha_upper_k - alpha_lower_k);
        alpha_middle_k = alpha_right_k;
        alpha_right_k = alpha_lower_k + (3/4) * (alpha_upper_k - alpha_lower_k);

        x_alpha_left_k = x_current + alpha_left_k * d_current;
        x_alpha_middle_k = x_current + alpha_middle_k * d_current;
        x_alpha_right_k = x_current + alpha_right_k * d_current;

        f_alpha_left_k = f_test(x_alpha_left_k);
        f_alpha_middle_k = f_alpha_right_k;
        f_alpha_right_k = f_test(x_alpha_right_k);
    end
k = k + 1;
end

alpha_star = (alpha_left_k + alpha_right_k) / 2;
x_next = x_current + alpha_star * d_current;
f_next = f_test(x_next);
end

2.主函数运行文件

放开相应注释即可。

% 这个是三点等间隔搜索法的主程序

% 清空变量
close;
clear;
clc;

% 示例一
% x_current = -0.5;
% d_current = 1;
% alpha_lower = 0;
% alpha_upper = 1;
% tolerance = 1e-4;
% [alpha_star, x_next, f_next, k] = Trisection_search(@f_test1, x_current, d_current, alpha_lower, alpha_upper, tolerance);

% 示例二
% x_current = [2, 2];
% d_current = [-1, -1];
% alpha_lower = 0;
% alpha_upper = 2;
% tolerance = 1e-6;
% [alpha_star, x_next, f_next, k] = Trisection_search(@f_test2, x_current, d_current, alpha_lower, alpha_upper, tolerance);

% 示例三
x_current = 3;
d_current = -1;
alpha_lower = 0;
alpha_upper = 2;
tolerance = 1e-4;
[alpha_star, x_next, f_next, k] = Trisection_search(@f_test3, x_current, d_current, alpha_lower, alpha_upper, tolerance);

3.一些说明

在三点搜索法的实现过程中，目前放开注释的是比较符合书中思路的，在前面有一部分注释掉的实现方式，如图所示。

使用这一部分代码时候，也是可以运行成功的，但是计算效率不高，因为三点搜索法的一些条件没有利用到，重新计算函数值这些会浪费一定时间，但是代码简单，比较好理解，也好书写。

（四） Fibonacci搜索法

1.Fibonacci搜索法函数文件

function [alpha_star, x_next, f_next, k] = Fibonacci_search(f_test, x_current, d_current, alpha_lower, alpha_upper, tolerance)
% 该函数针对Fibonacci搜索法的实现
% 输入参数说明
% f_test, 目标函数
% x_current, x在向量空间中的当前点（已确定）
% d_current, f_test在向量空间中的搜索方向（已确定）
% alpha_lower, 从x_current出发，沿着d_current得到一个步长的单谷区间，该单谷区间的初始下届
% alpha_upper, 从x_current出发，沿着d_current得到一个步长的单谷区间，该单谷区间的初始上届
% tolerance, 最终区间的要求宽度，精度不能小于扰动量
    
% 输出参数说明
% alpha_star, 完成对分搜索后输出的步长
% x_next, x_next = x_current + alpha_star * d_current;局部最优点
% f_next, 局部最优点对应的函数值
% k,迭代次数

% 计算斐波那契数列数列需要多少个
Fibonacci_series_upper = (alpha_upper - alpha_lower) / tolerance;
% 这里不需要第1个，实际上是从第二个开始
Fibonacci_series = [1, 2];
n = 2;
while (Fibonacci_series(n) <= Fibonacci_series_upper)
    n = n + 1;
    Fibonacci_series(n) = Fibonacci_series(n-1) + Fibonacci_series(n-2);
end

% 使用斐波那契数列进行搜索，注意到第n-2位置
k = 0;
alpha_lower_k = alpha_lower;
alpha_upper_k = alpha_upper;
Length_k = (Fibonacci_series(n - 1) / Fibonacci_series(n)) * (alpha_upper_k - alpha_lower_k);
alpha_left_k = alpha_upper_k - Length_k;
alpha_right_k = alpha_lower_k + Length_k;
x_alpha_left_k = x_current + alpha_left_k * d_current;
x_alpha_right_k = x_current + alpha_right_k * d_current;
f_alpha_left_k = f_test(x_alpha_left_k);
f_alpha_right_k = f_test(x_alpha_right_k);
while k < n - 2
    k = k + 1;
    if(f_alpha_left_k < f_alpha_right_k)
        % alpha_lower_k = alpha_lower;
        alpha_upper_k = alpha_right_k;
        Length_k = (Fibonacci_series(n - k - 1) / Fibonacci_series(n - k)) * (alpha_upper_k - alpha_lower_k);
        alpha_right_k = alpha_left_k;
        alpha_left_k = alpha_upper_k - Length_k;
        % alpha_right_k = alpha_lower_k + Length_k;  
        x_alpha_left_k = x_current + alpha_left_k * d_current;
        x_alpha_right_k = x_current + alpha_right_k * d_current;
        f_alpha_left_k = f_test(x_alpha_left_k);
        f_alpha_right_k = f_test(x_alpha_right_k);
    else
        alpha_lower_k = alpha_left_k;
        % alpha_upper_k = alpha_right_k;
        Length_k = (Fibonacci_series(n - k - 1) / Fibonacci_series(n - k)) * (alpha_upper_k - alpha_lower_k);
        alpha_left_k = alpha_right_k;
        % alpha_left_k = alpha_upper_k - Length_k;
        alpha_right_k = alpha_lower_k + Length_k;  
        x_alpha_left_k = x_current + alpha_left_k * d_current;
        x_alpha_right_k = x_current + alpha_right_k * d_current;
        f_alpha_left_k = f_test(x_alpha_left_k);
        f_alpha_right_k = f_test(x_alpha_right_k);
    end
end

k = k + 1;
disturbance_quantity = 0.1 * tolerance;
alpha_left_k = (1/2)*(alpha_lower_k + alpha_upper_k) - disturbance_quantity;
alpha_right_k = (1/2)*(alpha_lower_k + alpha_upper_k) + disturbance_quantity;
x_alpha_left_k = x_current + alpha_left_k * d_current;
x_alpha_right_k = x_current + alpha_right_k * d_current;
f_alpha_left_k = f_test(x_alpha_left_k);
f_alpha_right_k = f_test(x_alpha_right_k);

if (f_alpha_left_k < f_alpha_right_k)
    alpha_upper_k = alpha_right_k;
elseif (f_alpha_left_k > f_alpha_right_k)
    alpha_lower_k = alpha_left_k;
else
    alpha_upper_k = alpha_right_k;
    alpha_lower_k = alpha_left_k;
end
 
alpha_star = (alpha_lower_k + alpha_upper_k) / 2;
x_next = x_current + alpha_star * d_current;
f_next = f_test(x_next);
end

2.主函数运行文件

放开相应注释即可。

% 这个文件主要用来进行斐波那契搜索法的运行

close;
clear;
clc;

% 示例一
% x_current = -0.5;
% d_current = 1;
% alpha_lower = 0;
% alpha_upper = 1;
% tolerance = 1e-4;
% [alpha_star, x_next, f_next, k] = Fibonacci_search(@f_test1, x_current, d_current, alpha_lower, alpha_upper, tolerance);

% 示例二
x_current = [2, 2];
d_current = [-1, -1];
alpha_lower = 0;
alpha_upper = 2;
tolerance = 1e-6;
[alpha_star, x_next, f_next, k] = Fibonacci_search(@f_test2, x_current, d_current, alpha_lower, alpha_upper, tolerance);

% 示例三
% x_current = 3;
% d_current = -1;
% alpha_lower = 0;
% alpha_upper = 2;
% tolerance = 1e-4;
% [alpha_star, x_next, f_next, k] = Fibonacci_search(@f_test3, x_current, d_current, alpha_lower, alpha_upper, tolerance);

（五）黄金分割法

1.黄金分割法函数文件

function [alpha_star, x_next, f_next, k] = Golden_section_search(f_test, x_current, d_current, alpha_lower, alpha_upper, tolerance)
% 该函数针对黄金分割法的实现
% 输入参数说明
% f_test, 目标函数
% x_current, x在向量空间中的当前点（已确定）
% d_current, f_test在向量空间中的搜索方向（已确定）
% alpha_lower, 从x_current出发，沿着d_current得到一个步长的单谷区间，该单谷区间的初始下届
% alpha_upper, 从x_current出发，沿着d_current得到一个步长的单谷区间，该单谷区间的初始上届
% tolerance, 最终区间的要求宽度，精度不能小于扰动量
    
% 输出参数说明
% alpha_star, 完成对分搜索后输出的步长
% x_next, x_next = x_current + alpha_star * d_current;局部最优点
% f_next, 局部最优点对应的函数值
% k,迭代次数

% 黄金分割点
golden_ratio = 1 / (0.5 * (1 + sqrt(5)));

% 使用斐波那契数列进行搜索，注意到第n-2位置
k = 0;
alpha_lower_k = alpha_lower;
alpha_upper_k = alpha_upper;
Length_k = golden_ratio * (alpha_upper_k - alpha_lower_k);
alpha_left_k = alpha_upper_k - Length_k;
alpha_right_k = alpha_lower_k + Length_k;
x_alpha_left_k = x_current + alpha_left_k * d_current;
x_alpha_right_k = x_current + alpha_right_k * d_current;
f_alpha_left_k = f_test(x_alpha_left_k);
f_alpha_right_k = f_test(x_alpha_right_k);
while abs(alpha_upper_k - alpha_lower_k) > tolerance
    k = k + 1;
    if(f_alpha_left_k < f_alpha_right_k)
        % alpha_lower_k = alpha_lower;
        alpha_upper_k = alpha_right_k;
        Length_k = golden_ratio * (alpha_upper_k - alpha_lower_k);
        alpha_right_k = alpha_left_k;
        alpha_left_k = alpha_upper_k - Length_k;
        % alpha_right_k = alpha_lower_k + Length_k;  
        x_alpha_left_k = x_current + alpha_left_k * d_current;
        x_alpha_right_k = x_current + alpha_right_k * d_current;
        f_alpha_left_k = f_test(x_alpha_left_k);
        f_alpha_right_k = f_test(x_alpha_right_k);
    else
        alpha_lower_k = alpha_left_k;
        % alpha_upper_k = alpha_right_k;
        Length_k = golden_ratio * (alpha_upper_k - alpha_lower_k);
        alpha_left_k = alpha_right_k;
        % alpha_left_k = alpha_upper_k - Length_k;
        alpha_right_k = alpha_lower_k + Length_k;  
        x_alpha_left_k = x_current + alpha_left_k * d_current;
        x_alpha_right_k = x_current + alpha_right_k * d_current;
        f_alpha_left_k = f_test(x_alpha_left_k);
        f_alpha_right_k = f_test(x_alpha_right_k);
    end
end
 
alpha_star = (alpha_lower_k + alpha_upper_k) / 2;
x_next = x_current + alpha_star * d_current;
f_next = f_test(x_next);
end

2.主函数运行文件

% 这个文件主要用来进行黄金分割法的运行

close;
clear;
clc;

% 示例一
% x_current = -0.5;
% d_current = 1;
% alpha_lower = 0;
% alpha_upper = 1;
% tolerance = 1e-4;
% [alpha_star, x_next, f_next, k] = Golden_section_search(@f_test1, x_current, d_current, alpha_lower, alpha_upper, tolerance);

% 示例二
x_current = [2, 2];
d_current = [-1, -1];
alpha_lower = 0;
alpha_upper = 2;
tolerance = 1e-6;
[alpha_star, x_next, f_next, k] = Golden_section_search(@f_test2, x_current, d_current, alpha_lower, alpha_upper, tolerance);

% 示例三
% x_current = 3;
% d_current = -1;
% alpha_lower = 0;
% alpha_upper = 2;
% tolerance = 1e-4;
% [alpha_star, x_next, f_next, k] = Golden_section_search(@f_test3, x_current, d_current, alpha_lower, alpha_upper, tolerance);

3.一些说明

这里我稍微更改了一下实现方式，主要是针对于更新后的边界。书中原有方法如下所示，由于黄金分割法跟Fibonacci搜索法基本原理相似，所以，我就采用跟Fibonacci搜索法相似的更新方法。

更新方法如下所示。写的比较粗糙，但原理一致。

（六）三点二次插值法

1.三点二次插值法

function [alpha_star, x_next, f_next, k] = Quadratic3points_search(f_test, x_current, d_current, alpha_lower, alpha_upper, tolerance)
% 该函数针对三点二次插值法的实现
% 输入参数说明
% f_test, 目标函数
% x_current, x在向量空间中的当前点（已确定）
% d_current, f_test在向量空间中的搜索方向（已确定）
% alpha_lower, 从x_current出发，沿着d_current得到一个步长的单谷区间，该单谷区间的初始下届
% alpha_upper, 从x_current出发，沿着d_current得到一个步长的单谷区间，该单谷区间的初始上届
% tolerance, 最终区间的要求宽度，精度不能小于扰动量
    
% 输出参数说明
% alpha_star, 完成对分搜索后输出的步长
% x_next, x_next = x_current + alpha_star * d_current;局部最优点
% f_next, 局部最优点对应的函数值
% k,迭代次数

 
    
% 初始值
k = 0;
alpha_lower_k = alpha_lower;
alpha_upper_k = alpha_upper;

alpha_left_k = alpha_lower_k;
alpha_right_k = alpha_upper_k;
alpha_middle_k = (alpha_left_k + alpha_right_k) / 2;

x_alpha_left_k = x_current + alpha_left_k * d_current;
x_alpha_right_k = x_current + alpha_right_k * d_current;
x_alpha_middle_k = x_current + alpha_middle_k * d_current;

f_alpha_left_k = f_test(x_alpha_left_k);
f_alpha_right_k = f_test(x_alpha_right_k);
f_alpha_middle_k = f_test(x_alpha_middle_k);

% 用来保证算法能够启动，因为需要|ak+1-ak|作为判断条件，但第一次无法
% 计算这个，通过100 * tolerance保证能够进入循环进行运算，之后再进行
% 更新
gap_k = 100 * tolerance;
alpha_interpolation_previous = -100;

while gap_k > tolerance
    k = k + 1;
    a1 = (f_alpha_left_k - f_alpha_middle_k) * (alpha_middle_k - alpha_right_k) * (alpha_right_k - alpha_left_k);
    a2 = f_alpha_left_k * (alpha_middle_k - alpha_right_k) + f_alpha_middle_k * (alpha_right_k - alpha_left_k) + f_alpha_right_k * (alpha_left_k - alpha_middle_k);
    alpha_interpolation_k = (1/2) * ((alpha_left_k + alpha_middle_k) + (a1 / a2));
    x_alpha_interpolation_k = x_current + alpha_interpolation_k * d_current;
    f_alpha_interpolation_k = f_test(x_alpha_interpolation_k);

    alpha_interpolation_current = alpha_interpolation_k;
    alpha_middle_previous = alpha_middle_k;
    f_alpha_middle_previous = f_alpha_middle_k;

    if ((alpha_left_k < alpha_interpolation_k) && (alpha_interpolation_k < alpha_middle_k))
        if (f_alpha_interpolation_k <= f_alpha_middle_k)
            alpha_right_k = alpha_middle_k;
            f_alpha_right_k = f_alpha_middle_k;
            alpha_middle_k = alpha_interpolation_k;
            f_alpha_middle_k = f_alpha_interpolation_k;
        else
            alpha_left_k = alpha_interpolation_k;
            f_alpha_left_k = f_alpha_interpolation_k;
        end
    end

    if ((alpha_middle_k < alpha_interpolation_k) && (alpha_interpolation_k < alpha_right_k))
        if (f_alpha_interpolation_k <= f_alpha_middle_k)
            alpha_left_k = alpha_middle_k;
            f_alpha_left_k = f_alpha_middle_k;
            alpha_middle_k = alpha_interpolation_k;
            f_alpha_middle_k = f_alpha_interpolation_k;
        else
            alpha_right_k = alpha_interpolation_k;
            f_alpha_right_k = f_alpha_interpolation_k;
        end
    end
    gap_k = abs(alpha_interpolation_current - alpha_interpolation_previous);
    alpha_interpolation_previous = alpha_interpolation_current;
end
    
f_alpha_interpolation_current = f_alpha_interpolation_k;
if (f_alpha_interpolation_current <= f_alpha_interpolation_k)
    alpha_star = alpha_interpolation_k;
else
    alpha_star = alpha_middle_previous;
end

x_next = x_current + alpha_star * d_current;
f_next = f_test(x_next);
end

2.主函数运行文件

% 这个文件主要用来进行斐波那契搜索法的运行

close;
clear;
clc;

% 示例一
% x_current = -0.5;
% d_current = 1;
% alpha_lower = 0;
% alpha_upper = 1;
% tolerance = 1e-4;
% [alpha_star, x_next, f_next, k] = Quadratic3points_search(@f_test1, x_current, d_current, alpha_lower, alpha_upper, tolerance);

% 示例二
% x_current = [2, 2];
% d_current = [-1, -1];
% alpha_lower = 0;
% alpha_upper = 2;
% tolerance = 1e-6;
% [alpha_star, x_next, f_next, k] = Quadratic3points_search(@f_test2, x_current, d_current, alpha_lower, alpha_upper, tolerance);

% 示例三
x_current = 3;
d_current = -1;
alpha_lower = 0;
alpha_upper = 2;
tolerance = 1e-4;
[alpha_star, x_next, f_next, k] = Quadratic3points_search(@f_test3, x_current, d_current, alpha_lower, alpha_upper, tolerance);

3.一些说明

如果按照文中代码，如图所示红色部分圈出来的。在针对第一个和第二个示例的时，运行结果跟书中一致，跟前面几种方法是相同的。

但如果运行第三个示例的话，按照书中这种方法，则不会运行处正确结果。原因是因为这样计算出来的gap_k永远为0，程序只会迭代一次就结束。

更改完成的代码如上面所示，主要是提前将alpha_interpolation_previous进行初始化，这个值我设置为了-100，可以设置不可能达到的数字，这样能够保证每次运行都可以成功。其次，我将alpha_interpolation_previous的更新放在循环的最后，这样解决了上述问题。

GFCGUO

关注

8
点赞
踩
5

收藏

觉得还不错? 一键收藏
打赏
0
评论
应用最优化方法及MATLAB实现——第3章代码实现

在阅读最优方法及MATLAB实现后，想着将书中提供的代码自己手敲一遍，来提高自己对书中内容理解程度，巩固一下。这部分内容主要针对第3章的内容，将其所有代码实现均手敲一遍，中间部分代码自己根据其公式有些许的调整，但最后运算结果与书中所给结果一致。修改的部分我会在下面进行单独说明并标注上我这样做的原因。
复制链接

扫一扫