常用分式规划变换简述与仿真实验
在做课题时遇到一个子问题为线性分式规划问题，这里尝试用不同方法求解，同时做一些记录！
1. 线性分式规划问题模型
考虑的线性分式规划(Liner )问题定义如下：
\[\begin{align*} (\mathrm{P1})\quad & \\ \underset{\bm{\tau}}{\max}\quad &\frac{\bm{r}^T\bm{\tau}}{\bm{p}^T\bm{\tau}+p_c} \\ \mathrm{s.t.}\quad &\bm{\tau} \ge \bm{c}_1 \\ & \bm{C}_2\bm{\tau} \ge \bm{0} \\ & \bm{r}_d^T\bm{\tau} \ge c_3 \\ & \bm{1}^T\bm{\tau} \le 1 \end{align*} \]问题\((\mathrm{P1})\)的决策变量为 \(\bm{\tau} \in \mathbb{R}_+^K\)。其余各项参数有如下限制：\(\bm{r} > \bm{0}，\bm{p} > \bm{0}，p_c > 0，\bm{c}_1 > \bm{0}\)。此外，其中\(\bm{1}\)表示合适长度的全1向量。
2. Charnes-Cooper变换方法
Charnes-Cooper变换方法1 的核心思想有如下变换：
\[\begin{align*} z = {}& \frac{1}{\bm{p}^T\bm{\tau}+p_c} \\ \bm{q} ={}& \frac{\bm{\tau}}{\bm{p}^T\bm{\tau}+p_c} \end{align*} \]其中: \(z\in \mathbb{R}_+，\bm{q} \in \mathbb{R}_+^K\)。
将原问题\((\mathrm{P1})\)中的决策变量\(\bm{\tau}\)用 \(\frac{\bm{q}}{z}\) 代替可转化为等价问题\((\mathrm{P2})\)：
\[\begin{align*} (\mathrm{P2})\quad & \\ \underset{z,\bm{q}}{\max}\quad &\bm{r}^T\bm{q} \\ \mathrm{s.t.}\quad &\bm{q} \ge z\bm{c}_1 \\ & \bm{C}_2\bm{q} \ge \bm{0} \\ & \bm{r}_d^T\bm{q} \ge zc_3 \\ & \bm{1}^T\bm{q} \le z \\ & \bm{p}^T\bm{q} + zp_c = 1 \\ & \bm{q} > \bm{0}，z >0 \end{align*} \]显然问题\((\mathrm{P2})\)是关于决策变量\(z,\bm{q}\)的线性规划问题，可以直接使用求解器高效求解求解。获得最优解\(z^*,\bm{q}^*\) 后，则原问题\((\mathrm{P1})\)的最优优解为：
\[\bm{\tau}^* = \frac{\bm{q}^*}{z^*} \]3. Dinkelbach变换方法
3.1. Dinkelbach变换方法介绍
Dinkelbach转换方法2 解决如下形式的分式规划问题：
\[\begin{align*} \underset{\bm{x}}{\max}\quad & \frac{A(x)}{B(x)} \\ \mathrm{s.t.}\quad &\bm{x} \in \mathcal{X} \end{align*} \]其中\(\mathcal{X}\)为\(n\)维度欧几里得空间 \(\mathbb{R}^n\) 的一个紧闭子集，且要求：
对 \(\forall \bm{x} \in \mathcal{X}\) 都有 \(A(x) \ge 0, B(x) > 0\) 成立。
根据Dinkelbach的理论：
若
当且仅当：
\[\max \{ A(x)-Q^*B(x)| \bm{x} \in \mathcal{X} \} = 0 \]3.2. Dinkelbach算法总结
据此定义如下问题\((\mathrm{P-Dink})\)：
\[\begin{align*} (\mathrm{P-Dink})\quad & \\ \underset{\bm{x}}{\max}\quad & A(x)-QB(x) \\ \mathrm{s.t.}\quad &\bm{x} \in \mathcal{X} \end{align*} \]为了便于叙述，记函数\(F(x,Q) = A(x) - QB(x)\)。
Dinkelbach算法归纳如下：
1. 令 Q = 0, d = 0.0001, iter = 0, iterMax = 50
2. while iter < 50
3. 求解问题(P-Dink)得到最优解 x', 最优目标函数值 F'
4. if F' <= d
5 x_opt = x'
6. break
7. end
8. 令 Q = A(x')/B(x')
9. iter = iter + 1
10. end
11. 输出结果: x_opt
3.3. 本例的Dinkelbach变换
本例的\(\mathrm{(P-Dink)}\)问题如下所述：
\[\begin{align*} (\mathrm{P-Dink-P1})\quad & \\ \underset{\bm{\tau}}{\max}\quad &\bm{r}^T\bm{\tau}-Q(\bm{p}^T\bm{\tau}+p_c) \\ \mathrm{s.t.}\quad &\bm{\tau} \ge \bm{c}_1 \\ & \bm{C}_2\bm{\tau} \ge \bm{0} \\ & \bm{r}_d^T\bm{\tau} \ge c_3 \\ & \bm{1}^T\bm{\tau} \le 1 \end{align*} \]易知问题\((\mathrm{P-Dink-P1})\)是线性规划问题，所以可以获得全局最优解；故而通过Dinkelbach算法交替迭代更新\(\bm{\tau},Q\)，最终会分式规划问题的收敛到全局最优解。
4. quadratic变换
4.1. quadratic变换介绍
quadratic变换3 主要用来求解多项比分式规划(multiple-ratio FP)问题的。这里就简单介绍参考文献 [3] 提出的多项比和最大化问题的二次变换。
定义如下多项比和最大化问题 \((\mathrm{MRS-FP})\)：
\[\begin{align*} (\mathrm{MRS-FP})\quad & \\ \underset{\bm{x}}{\max}\quad & \sum_{m=1}^M \frac{A_m(x)}{B_m(x)} \\ \mathrm{s.t.}\quad &\bm{x} \in \mathcal{X} \end{align*} \]其中 \(n\)维欧几里得子空间 \(\mathcal{X}\) 为决策变量 \(x\) 的可行域；并且要求对任意\(m\)，函数\(A_m(x)\)值非负数，函数\(A_m(x)\)为值严格正数，即\(A_m(x): \mathbb{R}^n\rightarrow \mathbb{R}_{+}，B_m(x): \mathbb{R}^n\rightarrow \mathbb{R}_{++}\)。
使用二次转换将问题 \((\mathrm{MRS-FP})\) 转换为问题 \((\mathrm{MRS-FP-Quad})\):
\[\begin{align*} (\mathrm{MRS-FP-Quad})\quad & \\ \underset{\bm{x},\bm{y}}{\max}\quad & \sum_{m=1}^M (2y_m\sqrt{A_m(x)} -y_m^2B_m(x)) \\ \mathrm{s.t.}\quad &\bm{x} \in \mathcal{X} \\ & \bm{y} \in \mathbb{R}^M \end{align*} \]在获得问题 \((\mathrm{MRS-FP})\) 后，与类似与Dinkelbach变换，quadratic变换也是交替更新变量 \(\bm{x}，\bm{y}\)，直到收敛。
4.2. Quadratic变换算法总结
求解上述问题的Quadratic变换算法总结如下图所示：
对于\(A_m(x)\)为凹函数，\(B_m(x)\)为凸函数，可行域\(\mathcal{X}\) 为凸集的情况，上述算法可以保障收敛最终收敛点为驻点；更进一步如果 \(M=1\) 即问题退化为 单项比分式规划问题(前面两种变换介绍的例子)，那么该算法可以收敛到全局最优解3。另外文献 [3] 还讨论了多项比函数和问题，多项比最小最大问题以及多维复空间情形的Quadratic变换形式。
4.3. 本例的Quadratic变换形式
对比发现本例问题\((\mathrm{P1})\)的 \(M=1\), \(A(\bm{\tau}) = \bm{r}^T\bm{\tau}\)，\(B(\bm{\tau}) = \bm{p}^T\bm{\tau} + p_c\)
那么其Quadratic变换形式为：
显而易见，当固定 \(y\) 时，问题\((\mathrm{Quadratic-P1})\) 是关于 \(\bm{\tau}\) 的凸优化问题，故而上述Quadratic变换算法的第4步可求得全局最优解；并且对于该问题，Quadratic变换算法可以收敛到全局最优解。
5. 数值计算及算法实现
5.1. 参数设定
给定参数如 ParamSet.m 脚本所下示：
% ParamSet.m脚本定义优化问题中的参数
%决策变量维度K
K = 3;
r =[20.8204;
24.8497;
22.5085];
p = [ 0.97;
1;
0.99];
p_c = 1;
c_1 = [ 0.0451;
0.0408;
0.0312];
C_2 = [ -0.0010 0.0493 0.0202;
0.0313 -0.0010 0.0202;
0.0313 0.0493 -0.0010];
r_d = [ 5.7464
8.2832
9.4225];
c_3 = 8;
5.2. Charnes-Cooper变换算法实现
Charnes-Cooper变换算法实现代码为如下 LFP_Charnes_Cooper.m 脚本所示：
function [isFeasible,obj_opt,tau_opt] = LFP_Charnes_Cooper(K,r,p,p_c,c_1,C_2,r_d,c_3)
%LFP_CHARNES_COOPER 使用Charnes-Cooper变换求解线性分式规划问题
% 输入参数说明：
% K,r,p,p_c,c_1,C_2,r_d,c_3 优化问题参数
% 输出参数说明：
% isFeasible 是否可行
% obj_opt 最优目标函数值
% tau_opt 最优解
%% 调用CVX求解线性规划
cvx_begin quiet
variable q(K,1) nonnegative
variable z nonnegative
maximize (r'*q)
subject to
q >= z * c_1;
C_2 * q >= 0;
r_d'* q >= z * c_3;
sum(q) <= z;
p'* q + z * p_c == 1;
cvx_end
%% 结果返回
obj_opt = cvx_optval;
isFeasible = true;
tau_opt = q/z;
if strcmp(cvx_status,'Infeasible') || strcmp(cvx_status,'Unbounded') || isnan(cvx_optval) || isinf(cvx_optval)
isFeasible = false;
end
end
5.3. Dinkelbach变换算法实现
Dinkelbach变换算法实现代码为如下 LFP_Dinkelbach.m 脚本所示：
function [isFeasible,obj_opt,tau_opt,Q_log,F_log] = LFP_Dinkelbach(K,r,p,p_c,c_1,C_2,r_d,c_3,delta,IterMax)
%LFP_DINKELBACH 使用Dinkelbach变换求解线性分式规划问题
% 输入参数说明：
% K,r,p,p_c,c_1,C_2,r_d,c_3 优化问题参数
% delta 一个小正数 收敛门限
% IterMax 最大迭代次数
% 输出参数说明：
% isFeasible 是否可行
% obj_opt 最优目标函数值
% tau_opt 最优解
% Q_log 记录每次迭代的Q值
% F_log 记录每次迭代的F值
%% 初始化
Q_log = [];
F_log = [];
isFeasible = true;
Q = 0;
%% 迭代求解
for j = 1 : IterMax
cvx_clear
% 调用CVX
cvx_begin quiet
variable tau(K,1) nonnegative
maximize ( r'*tau - Q *(p'*tau + p_c) )
subject to
tau >= c_1;
C_2 * tau >= 0;
r_d'* tau >= c_3;
sum(tau) <= 1;
cvx_end
% 不可行或者无界则直接退出循环
if strcmp(cvx_status,'Infeasible') || strcmp(cvx_status,'Unbounded') || isnan(cvx_optval) || isinf(cvx_optval)
isFeasible = false;
break
end
Q_log = [Q_log;Q];
F_log = [F_log;cvx_optval];
% 如果收敛则退出循环，如果不收敛则更新Q
if cvx_optval <= delta
break
else
Q = r'*tau / (p'*tau + p_c);
end
end
%% 返回结果
tau_opt = tau;
obj_opt = r'*tau_opt / (p'*tau_opt + p_c);
end
5.4. Quadratic变换算法实现
Quadratic变换算法实现代码为如下 LFP_Quadratic.m 脚本所示：
function [isFeasible,obj_opt,tau_opt,F_log] = LFP_Quadratic(K,r,p,p_c,c_1,C_2,r_d,c_3,delta,IterMax)
%LFP_QUADRATIC 使用Quadratic变换求解线性分式规划问题
% 输入参数说明：
% K,r,p,p_c,c_1,C_2,r_d,c_3 优化问题参数
% delta 一个小正数 收敛门限
% IterMax 最大迭代次数
% 输出参数说明：
% isFeasible 是否可行
% obj_opt 最优目标函数值
% tau_opt 最优解
% F_log 记录每次迭代的F值
%% 计算初始可行解
% 使用CVX求初始解 一个线性规划可行性问题
cvx_begin quiet
variable tau(K,1) nonnegative
subject to
tau >= c_1;
C_2 * tau >= 0;
r_d'* tau >= c_3;
sum(tau) <= 1;
cvx_end
% 初始化数据
tau_i = tau;
F_i = r'*tau_i / (p'*tau_i + p_c);
F_log = [F_i];
isFeasible = true;
%% 迭代求解
%% 迭代求解
for j = 1 : IterMax
y = sqrt(r'*tau_i)/(p'*tau_i + p_c);
cvx_clear
% 调用CVX
cvx_begin quiet
variable tau(K,1) nonnegative
maximize ( 2*y*sqrt(r'*tau) - y^2*(p'*tau + p_c) )
subject to
tau >= c_1;
C_2 * tau >= 0;
r_d'* tau >= c_3;
sum(tau) <= 1;
cvx_end
tau_i_1 = tau;
F_i_1 = r'*tau_i_1 / (p'*tau_i_1 + p_c);
F_log = [F_log;F_i_1];
% 不可行或者无界则直接退出循环
if strcmp(cvx_status,'Infeasible') || strcmp(cvx_status,'Unbounded') || isnan(cvx_optval) || isinf(cvx_optval)
isFeasible = false;
break
end
% 如果收敛则退出循环，如果不收敛则更新Q
if (F_i_1 - F_i) <= delta
break
else
tau_i = tau_i_1;
F_i = F_i_1;
end
end
obj_opt = F_i_1;
tau_opt = tau_i_1;
end
5.5. 结果对比
由于本身这个线性分式规划就比较简单，维度有比较小，所有在用Dinkelbach变换算法和Quadratic变换算法迭代求解时，一两步就收敛了所以就不展示其收敛曲线了。
然后为了观察，对优化问题中的参数 \(c_3\) 进行参数扫描。仿真代码为 Test.m脚本如下所示：
% Test.m 脚本
clc
ParamSet
delta = 1E-6;
IterMax = 20;
%%
c_3_vec = 7.6:0.2:9.2;
L = length(c_3_vec);
result_CC = zeros(L,1);
result_Dinkelbach = zeros(L,1);
result_Quadratic = zeros(L,1);
for j = 1 : L
[isFeasible,obj_opt,tau_opt_CC] = LFP_Charnes_Cooper(K,r,p,p_c,c_1,C_2,r_d,c_3_vec(j));
result_CC(j) = obj_opt;
[isFeasible,obj_opt,tau_opt_Dink,Q_log,F_log] = LFP_Dinkelbach(K,r,p,p_c,c_1,C_2,r_d,c_3_vec(j),delta,IterMax);
result_Dinkelbach(j) = obj_opt;
[isFeasible,obj_opt,tau_opt_Quad,F_log] = LFP_Quadratic(K,r,p,p_c,c_1,C_2,r_d,c_3_vec(j),delta,IterMax);
result_Quadratic(j) = obj_opt;
fprintf('# %d \n %s \n %s \n %s \n',j, num2str(tau_opt_CC'), num2str(tau_opt_Dink'), num2str(tau_opt_Quad'))
end
%%
plot(c_3_vec,result_CC,'-o',c_3_vec,result_Dinkelbach,'-p',c_3_vec,result_Quadratic,'-+','MarkerSize',8,'LineWidth',2)
legend({'Charnes-Cooper','Dinkelbach','Quadratic'})
ylabel('最优目标值')
xlabel('c_3')
title('不同变换方法求解的最优目标函数值与参数c_3的关系')
绘制出各种变换算法的 最优目标值和参数 \(c_3\) 的变化曲线，如下图所示：
三条曲线完全重合，这是符号预期的，因为三种算法对于该线性分式规划问题都可以求解得全局最优解。
6. 总结
Charnes-Cooper变换可以将 单项比线性分式规划 (就是\((\mathrm{P1})\) 这种形式的问题) 转换为一个线性规划问题，不需要迭代，直接求解效率最高。如果用 Dinkelbach变换 和 Quadratic变换，就有点杀鸡用牛刀的意思了。但是对于分子分母非线性的情形，Charnes-Cooper变换引入了新的一个约束不一定好处理，并且可能变换后也不好求解，但具体问题具体分析，多试试。
Dinkelbach变换的适用性跟广泛，也不会引入新的约束。特别是如果分子项为凹函数，分母项为凸函数，可行域为凸集的情形，Dinkelbach变换后的问题为凸优化问题，可以高效求解，此时，通过迭代可收敛到全局最优解。
Quadratic变换的主要目的是简化 多项比线性分式规划 的求解，对于 单项比线性分式规划 的情形，在一些case下，收敛速度明显慢于Dinkelbach变换方法3。此外，考虑Quadratic变换后的目标函数含有 根号项，在一些问题中可能会变得不好处理。