手动求导真的很麻烦,针对之前上传的matlab代码,需要自己手算求梯度和黑森矩阵,所以改个代码不需要自己求,要matlab给我求。(所以说懒人推动科技发展bushi)
% % 定义目标函数 Beale 函数
% f = @(x)(x(1)-10^6)^2+(x(2)-2*10^-6)^2+(x(1)*x(2)-2)^2;
%
% % 初始点设定
% x0 = [1; 1];
%
% % 设置优化选项
% options = optimoptions(@fminunc,'Display','iter','Algorithm','quasi-newton','MaxIterations',3000);
% options = optimoptions(@fminunc,'Display','iter','Algorithm','quasi-newton','MaxIterations',3000, ...
% 'TolFun', 1e-6, 'TolX', 1e-6); % 设置解的变化的终止条件,当解向量的元素变化小于这个值时停止迭代
% % 使用fminunc函数求解最小值
% [x, fval, exitflag, output] = fminunc(f, x0, options);
%
% % 输出结果
% disp(['迭代次数: ', num2str(output.iterations)]);
% disp(['x的最优解: ', num2str(x')]);
% disp(['f的最小值: ', num2str(fval)]);
% % 绘制迭代图像
% figure;
% plot(1:output.iterations, output.fval, 'bo-');
% xlabel('Iterations');
% ylabel('Objective Function Value');
% title('Convergence of Objective Function Value');
% 定义目标函数 Beale 函数
% f = @(x)2+2i-exp[i*x(1)]-exp[i*x(2)];
%
% % 初始点设定
% x0 = [1; 1];
%
% % 设置优化选项,包括函数值和解的精度
% options = optimoptions(@fminunc,'Display','iter','Algorithm','quasi-newton','MaxIterations',3000, ...
% 'TolFun', 1e-6, 'TolX', 1e-6); % 设置解的变化的终止条件,当解向量的元素变化小于1e-6时停止迭代
%
% % 使用fminunc函数求解最小值
% [x, fval, exitflag, output] = fminunc(f, x0, options);
%
% 输出结果
% disp(['迭代次数: ', num2str(output.iterations)]);
% disp(['x的最优解: ', num2str(x')]);
% disp(['f的最小值: ', num2str(fval)]);
%定义目标函数 Beale 函数
% f = @(x) (1.5 - x(1) + x(1)*x(2))^2 + (2.25 - x(1) + x(1)*x(2)^2)^2 + (2.625 - x(1) + x(1)*x(2)^3)^2;
%
% % 初始点设定
% x0 = [1; 1];
%
% % 设置优化选项,包括函数值和解的精度
% options = optimoptions(@fminunc,'Display','iter','Algorithm','quasi-newton','MaxIterations',3000, ...
% 'TolFun', 1e-6, 'TolX', 1e-6); % 设置解的变化的终止条件,当解向量的元素变化小于1e-6时停止迭代
%
% % 使用fminunc函数求解最小值
% [x, fval, exitflag, output] = fminunc(f, x0, options);
%
% % 输出结果
% disp(['迭代次数: ', num2str(output.iterations)]);
% disp(['x的最优解: ', num2str(x')]);
% disp(['f的最小值: ', num2str(fval)]);
% 定义Jennrich-Sampson函数
f = @(x) sum((2 + 2*(1:10) - (exp((1:10) * x(1)) + exp((1:10) * x(2)))).^2);
% 初始点
x0 = [0.3, 0.4];
% 设置优化选项
options = optimset('Display', 'iter', 'TolFun', 1e-8, 'MaxIter', 1000);
% 调用fminsearch进行优化
[x, fval, exitflag, output] = fminsearch(f, x0, options);
disp(['找到的极小点坐标为:', mat2str(x), ', 对应的函数值为:', num2str(fval)]);
% f = @(x) ((10^4 * x(1) * x(2) - 1)^2 + (exp(-x(1)) + exp(-x(2)) - 1.0001))^2;
%
% % 初始点设定
% x0 = [0; 1];
%
% % 设置优化选项,包括函数值和解的精度
% options = optimoptions(@fminunc,'Display','iter','Algorithm','quasi-newton','MaxIterations',3000, ...
% 'TolFun', 1e-6, 'TolX', 1e-6); % 设置解的变化的终止条件,当解向量的元素变化小于1e-6时停止迭代
%
% % 使用fminunc函数求解最小值
% [x, fval, exitflag, output] = fminunc(f, x0, options);
%
% % 输出结果
% disp(['迭代次数: ', num2str(output.iterations)]);
% disp(['x的最优解: ', num2str(x')]);
% disp(['f的最小值: ', num2str(fval)]);
% % 定义Jennrich-Sampson函数
% f = @(x) (1000*(x(3)-5/pi*atan(x(2)/x(1))^3+1000*(x(1)^2+x(2)^2)^1/2-1)^3+x(3)^3);
%
% % 初始点
% x0 = [-1,0,0];
%
% % 设置优化选项
% options = optimset('Display', 'iter', 'TolFun', 1e-8, 'MaxIter', 1000);
%
% % 调用fminsearch进行优化
% [x, fval, exitflag, output] = fminsearch(f, x0, options);
%
% disp(['找到的极小点坐标为:', mat2str(x), ', 对应的函数值为:', num2str(fval)]);
% % 定义目标函数
% syms x1 x2 x3
% f = 1000*(x3-5/pi*atan(x2/x1))^3 + 1000*((x1^2+x2^2)^(1/2)-1)^3 + x3^3;
%
% % 计算梯度
% grad_f = gradient(f, [x1, x2, x3]);
%
% % 将梯度转换为一维向量形式
% grad_f_vec = reshape(grad_f, [], 1);
%
% % 初始化为数值形式以便后续迭代
% x0 = [-1; 0; 0]; % 修改为一维向量形式
%
% alpha = 0.01; % 步长
% max_iter = 1000; % 最大迭代次数
% tolerance = 1e-6; % 公差
%
% % 迭代过程
% for i = 1:max_iter
% % 计算当前点的符号梯度,并转化为数值梯度
% symbolic_gradient = subs(grad_f_vec, [x1; x2; x3], x0);
% numeric_gradient = double(symbolic_gradient);
%
% % 当 x1 < 0 时,计算修正后的目标函数值,但梯度保持不变
% if x0(1) < 0
% f_modified = 1000*(x3-5/pi*atan(x2/x1)+0.5)^3 + 1000*((x1^2+x2^2)^(1/2)-1)^3 + x3^3;
% current_f = double(subs(f_modified, [x1; x2; x3], x0));
% else
% current_f = double(subs(f, [x1; x2; x3], x0));
% end
%
% % 更新解
% x0 = x0 - alpha * numeric_gradient;
%
% % 检查梯度范数是否小于公差
% gradient_norm = norm(numeric_gradient);
% if gradient_norm < tolerance
% break;
% end
% end
%
% % 输出最优解
% disp(['最优解为:', num2str(x0')]);
% % 定义目标函数 Beale 函数
% f = @(x)(x(1)-10^6)^2+(x(2)-2*10^-6)^2+(x(1)*x(2)-2)^2;
%
% % 初始点设定
% x0 = [1; 1];
%
% % 设置优化选项
% options = optimoptions(@fminunc,'Display','iter','Algorithm','quasi-newton','MaxIterations',3000, ...
% 'TolFun', 1e-6, 'TolX', 1e-6); % 设置解的变化的终止条件,当解向量的元素变化小于这个值时停止迭代
%
% % 创建空数组用于存储每次迭代的目标函数值
% func_values = [];
%
% % 自定义输出函数,用于保存每次迭代的目标函数值
% options.OutputFcn = @(x, optimValues, state) saveFuncValues(x, optimValues, state);
%
% % 使用fminunc函数求解最小值
% [x, fval, exitflag, output] = fminunc(f, x0, options);
%
% % 输出结果
% disp(['迭代次数: ', num2str(output.iterations)]);
% disp(['x的最优解: ', num2str(x')]);
% disp(['f的最小值: ', num2str(fval)]);
%
% % 绘制迭代过程中连续下降搜索的图像
% figure;
% plot(1:length(func_values), func_values, 'bo-');
% xlabel('Iterations');
% ylabel('Objective Function Value');
% title('Convergence of Objective Function Value');
%
% % 保存每次迭代的目标函数值
% function stop = saveFuncValues(x, optimValues, state)
% persistent func_values
% if strcmp(state, 'init')
% func_values = [];
% end
% func_values = [func_values; optimValues.fval];
% stop = false;
% end
>> untitled3
Iteration Func-count min f(x) Procedure
0 1 4171.31
1 3 4171.31 initial simplex
2 5 2369.05 expand
3 7 1218.17 expand
4 9 574.632 expand
5 11 163.782 expand
6 12 163.782 reflect
7 13 163.782 reflect
8 14 163.782 reflect
9 16 149.863 reflect
10 18 149.863 contract inside
11 20 126.375 expand
12 22 126.375 contract inside
13 24 125.574 contract inside
14 26 125.574 contract inside
15 28 125.574 contract inside
16 29 125.574 reflect
17 31 124.709 contract inside
18 33 124.709 contract outside
19 35 124.588 contract inside
20 37 124.532 contract inside
21 39 124.441 contract inside
22 41 124.441 contract inside
23 43 124.382 reflect
24 45 124.382 contract inside
25 47 124.372 reflect
26 48 124.372 reflect
27 50 124.367 contract inside
28 52 124.363 contract inside
29 54 124.363 contract inside
30 56 124.363 contract outside
31 58 124.362 contract inside
32 60 124.362 contract inside
33 62 124.362 contract inside
34 64 124.362 contract outside
35 66 124.362 contract inside
36 68 124.362 contract inside
37 70 124.362 reflect
38 72 124.362 contract inside
39 74 124.362 contract outside
40 76 124.362 contract inside
41 78 124.362 contract inside
42 80 124.362 contract inside
43 82 124.362 contract inside
44 84 124.362 reflect
45 86 124.362 contract inside
46 88 124.362 contract inside
47 90 124.362 contract inside
优化已终止:
当前的 x 满足使用 1.000000e-04 的 OPTIONS.TolX 的终止条件,
F(X) 满足使用 1.000000e-08 的 OPTIONS.TolFun 的收敛条件
找到的极小点坐标为:[0.257826180021054 0.257823437205372], 对应的函数值为:124.3622
运行结果为Jennrich-Sampson函数的最优点和对应的函数值,是正确的。不过其中也有一些函数不适用,必须要手动求导,大家可以自己试一试。也可以套用函数模板,去掉百分号都可以运行出来!!!(终于又可以删除一个文件了!!!!!)最速下降法改编的代码比牛顿法改编的要适应性更强一些,结果输出的正确性也相对高一点,估计也是牛顿法要求太多了(尤其是正定性)。还有最后一段代码,我是想把迭代图表示出来的,可惜并没有(本学渣能力真的有限),希望大家可以能有所启发,顺便把代码发评论区,嘿嘿嘿。