在使用yalmip和cplex求解器求解问题中遇到了目标函数包含非线性的问题,由于是初学者没有使用过线性化方法,在这里写下笔记作为记录方便后续的学习。
一阶泰勒展开公式:
其中L(x)为函数f(x)在点[a,f(a)]处的一节泰勒展开,一阶泰勒展开公式可以在点[a,f(a)]附近近似非线性函数f(x)。
上面是泰勒展开的方法,将上述方法扩展至二元函数可得到一下公式:
即:
因此,获得泰勒展开式时可分为以下步骤:
1、求展开点函数值 2、求展开点梯度 3、构造线性函数
代码如下:
% 中心点 i
center = a1(i, :);
a1_i = center(1);
a2_i = center(2);
% 计算函数在中心点的值
f_a = subs(f, {x1, x2}, {a1_i, a2_i});
% 计算函数的梯度
grad_f = gradient(f, [x1, x2]);
% 计算梯度在中心点的值
grad_f_a = subs(grad_f, {x1, x2}, {a1_i, a2_i});
% 构造线性化函数
L{i} = f_a + grad_f_a(1)*(x1 - a1_i) + grad_f_a(2)*(x2 - a2_i);
在实际使用过程中由于我需要在较大范围内使用此线性化函数,由于当选取的点离展开点较远时误差会增大,因此将我的函数定义域划分为多区域,之后在每个区域内选取中心点作为展开点,最终得到一个分区域多段线性化函数,举例代码如下:
% 1. 定义二元函数
syms x1 x2
f = x1^2 + 2*x2^2;
% 2. 划分范围
% 子区域的中心点
a1 = [0.5, 0.5; 1.5, 0.5; 0.5, 1.5; 1.5, 1.5];
regions = {'区域1', '区域2', '区域3', '区域4'};
% 3. 在每个子区域内进行线性化
n_regions = size(a1, 1);
L = cell(n_regions, 1); % 存储每个子区域的线性化函数
for i = 1:n_regions
% 选择子区域的中心点
center = a1(i, :);
a1_i = center(1);
a2_i = center(2);
% 计算函数在中心点的值
f_a = subs(f, {x1, x2}, {a1_i, a2_i});
% 计算函数的梯度
grad_f = gradient(f, [x1, x2]);
% 计算梯度在中心点的值
grad_f_a = subs(grad_f, {x1, x2}, {a1_i, a2_i});
% 构造线性化函数
L{i} = f_a + grad_f_a(1)*(x1 - a1_i) + grad_f_a(2)*(x2 - a2_i);
end
为了让线性化结果可以直观表现,将原函数和线性化函数以图形形式展示,代码如下:
% 4. 显示结果
for i = 1:n_regions
fprintf('子区域 %d (%s):\n', i, regions{i});
disp(['线性化函数: ', char(L{i})]);
end
% 5. 分段线性化函数的可视化
% 生成网格数据
[x1_grid, x2_grid] = meshgrid(linspace(0, 2, 50), linspace(0, 2, 50));
L_values = zeros(size(x1_grid));
for i = 1:n_regions
% 计算当前区域的线性化函数值
center = a1(i, :);
a1_i = center(1);
a2_i = center(2);
% 计算当前区域的线性化函数值
L_i = double(subs(L{i}, {x1, x2}, {x1_grid, x2_grid}));
% 确定当前区域
mask = (x1_grid >= a1_i - 0.5 & x1_grid < a1_i + 0.5) & ...
(x2_grid >= a2_i - 0.5 & x2_grid < a2_i + 0.5);
% 将线性化函数值应用到相应区域
L_values(mask) = L_i(mask);
end
% 绘制图形
figure;
subplot(1, 2, 1);
surf(x1_grid, x2_grid, double(subs(f, {x1, x2}, {x1_grid, x2_grid})));
title('原函数 f(x_1, x_2)');
xlabel('x_1');
ylabel('x_2');
zlabel('f(x_1, x_2)');
subplot(1, 2, 2);
surf(x1_grid, x2_grid, L_values);
title('分段线性化函数');
xlabel('x_1');
ylabel('x_2');
zlabel('L(x_1, x_2)');
在实际使用中,将线性化函数的变量改为我的实际应用变量:
for i = 1:n_regions
disp('构造线性函数区域');
disp(i);
% 选择子区域的中心点
center = zxpoint(i, :);
a1_i = center(1);
a2_i = center(2);
% 计算函数在中心点的值
Dis_zx = double(subs(Dis, {X_1, X_2}, {a1_i, a2_i}));
% 计算函数的梯度
grad_Dis = gradient(Dis, [X_1, X_2]);
% 计算梯度在中心点的值
grad_f_Dis = double(subs(grad_Dis, {X_1, X_2}, {a1_i, a2_i}));
% 构造线性化函数
Line{i} = @(Bianliang1, Bianliang2) Dis_zx + grad_f_Dis(1) * (Bianliang1 - a1_i) + grad_f_Dis(2) * (Bianliang2 - a2_i);
end
% 定义目标函数
Dis_values = sdpvar(n_regions, 1);
for i = 1:n_regions
Dis_values(i) = Line{i}(Bianliang1, Bianliang2);
end
Bianliang1和Bianliang2为举例中我的实际应用变量,在使用过程中替换为自己真实变量即可。
本人没有系统学习过MATLAB、yalmip、cplex和线性化方法,以上内容是自己的想法,想法简单不够严谨,请各位大佬批评指正。