笔记：常见的约束问题求解算法——乘子法和Frank-Wolfe算法

前言
本文介绍了约束最优化问题中常用的两种算法，乘子法和 Frank-Wolfe 算法。
最后通过matlab编程实现了以上两种算法，并对实际问题进行求解。

符号定义
假设约束最优化问题
min⁡f(x),其中(x∈Rn)\min f(x), 其中(x \in R^n)minf(x),其中(x∈Rn)，
s.t.gi(x)≥0,(i∈I={1,2,⋯,m1})s.t.\ g_i(x) \geq 0, (i \in I=\{1,2,\cdots, m_1\})s.t. gi(x)≥0,(i∈I={1,2,⋯,m1})
hj(x)=0,(j∈E={m1+1,m1+2,⋯,m}))h_j(x) = 0, (j \in E = \{m_1+1, m_1+2, \cdots, m\}))hj(x)=0,(j∈E={m1+1,m1+2,⋯,m}))
记函数f(x)的一阶导数为∇f(x)\nabla f(x)∇f(x)，二阶导数为∇2f(x)\nabla^2 f(x)∇2f(x).

一、乘子法

简介
乘子法(multiplier method)是一种约束极小化的算法，通过引入Lagrange函数来求解得到最优解的。
适用范围：无约束最优化问题
算法步骤
步0：取初始点x(0)∈Rnx^{(0)} \in R^nx(0)∈Rn，初始乘子向量λ(0)\lambda^{(0)}λ(0). 给定罚参数序列{μk}\{\mu_k \}{μk}，精度ϵ>0\epsilon>0ϵ>0. 令k=0.
步1：下式构造增广Lagrange函数Lμ(x,λ)L_{\mu} (x, \lambda)Lμ(x,λ).
Lμ(x,λ)=f(x)+12μ−1∑i∈I[min2{μgi(x)−λi,0}−λi2]−∑j∈Eλihi(x)+12μ∑j∈Ehi2(x)L_{\mu}(x, \lambda) = f(x) + \frac{1}{2}\mu^{-1}\sum_{i \in I}[min^2 \{\mu g_i(x)-\lambda_i, 0\}-\lambda_i^2] - \sum_{j \in E}\lambda_ih_i(x) + \frac{1}{2}\mu\sum_{j \in E} h_i^2(x)Lμ(x,λ)=f(x)+21μ−1∑i∈I[min2{μgi(x)−λi,0}−λi2]−∑j∈Eλihi(x)+21μ∑j∈Ehi2(x)
步2：以x(k−1)x^{(k-1)}x(k−1)作为初始点(k=0时，初始点任意)，求解无约束问题min⁡Lμk(x,λ(k))\min L_{\mu_k}(x, \lambda^{(k)})minLμk(x,λ(k))得解x(k)x^{(k)}x(k).
步3：若∣∣hE(x(k))∣∣+∣∣min⁡{gI(x(k)),μk−1λI(k)}∣∣≤ϵ||h_E(x^{(k)})|| + ||\min \{g_I(x^{(k)}), \mu^{-1}_k\lambda_I^{(k)}\}|| \leq \epsilon∣∣hE(x(k))∣∣+∣∣min{gI(x(k)),μk−1λI(k)}∣∣≤ϵ则得解x(k)x^{(k)}x(k).
步4：在下式中取x=x(k),λ=λ(k)x = x^{(k)}, \lambda = \lambda^{(k)}x=x(k),λ=λ(k)得λ(k+1)\lambda^{(k+1)}λ(k+1)，令k=k+1, 转步1.
λj+={λj−μhj(x),(j∈E)max⁡{λi−μigi(x),0},(i∈I)\lambda_j^+ = \left\{ \begin{array}{ll} %该矩阵一共3列，每一列都居中放置 \lambda_j - \mu h_j(x), & (j \in E)\\ \max \{\lambda_i - \mu_ig_i(x), 0\}, & (i \in I)\\ \end{array} \right. λj+={λj−μhj(x),max{λi−μigi(x),0},(j∈E)(i∈I)
其中，x,λ,μx, \lambda, \mux,λ,μ表示当前迭代点的值，λj+\lambda^+_jλj+表示下一次迭代的乘子向量.
matlab实现
说明：函数中使用的initpt函数（用于初始化x0x_0x0）是为了自己编写的辅助函数，这几个函数的具体实现可见最下方的链接处。

function [iter_num, time, epsilonk] = method_of_multipliers(nprob, n, x, epsilon, max_iter)
% 乘子法   2020.06.25
% @author: 豆奶
% 函数功能：实现乘子法     % 比较两个字符串 strcmp(a, 'De_Jong')
% 输入：
% nprob: 问题编号，
% n: 维度
% x: 自变量x
% epsilon: int,
% max_iter: 最大迭代次数
% 输出：
% iter_num: 迭代次数
% time: 计算时间
% epsilonk: KKT系统模的最终值t1=clock;  % 计时开始
% 步0：取初始点x，初始乘子向量. 给定罚参数序列，精度
if nargin==2x = initpt(n, nprob);epsilon = 1e-5;max_iter = 100;
end
mu = 0.5;
xk = x;
% 设置目标函数，约束函数
if strcmp(nprob, 'De_Jong')f = @(x) sum(x.^2);g1 =@(x) x+5.12;g2 =@(x) 5.12-x;g =@(x) [g1(x); g2(x)];h =@(x) 0;lambda_g = 10*ones(2*n, 1);  % (2n, 1)   和约束条件的个数有关lambda_h = 0;
elseif strcmp(nprob, 'Rastrigin')f = @(x)10*n + sum(x.^2 - 10*cos(2*pi*x));  g1 =@(x) x+5.12;g2 =@(x) 5.12-x;g =@(x) [g1(x); g2(x)];h =@(x) 0;lambda_g = ones(2*n, 1);  % (2n, 1)   和约束条件的个数有关lambda_h = 0;
end
for k=1:max_iter % 步1：由（12.27）构造增广Lagrange函数L.L =@(x) f(x) + 0.5.*(1/mu)*sum(min(mu*g(x)-lambda_g, 0).^2-lambda_g.^2) - sum(lambda_h'*h(x)) + 0.5*mu*sum(h(x).^2);% 步2：以x_{k-1}作为初始点（k=0时，初始点任意），求解无约束问题L得解x_k.xk = fmincon(L, xk);% 步3：判断是否小于epsilon，以跳出循环.eps_k = norm(h(xk)) + norm(min(g(xk), 1/mu*lambda_g));if eps_k<=epsilonbreak;end % 步4：在（12.28）中取x=xk, lambda=lambda_klambda_g = max(lambda_g-mu*g(xk), 0);  % update lambda_g, i \in Ilambda_h = lambda_h - mu*h(xk);        % update lambda_h, j \in E
end
% 获取输出值
iter_num = k;
syms x;                     %通过符号变量将匿名函数转换为符号函数
gxk = mu*g(xk) - lambda_g;  % 由于min无法直接求导，因此将这部分割裂开来，独自求导计算
%dg = diff(g(x));  % 由于本次测试采用的两个测试问题的g的导数均为实数，即(1,-1)，失去了自变量x，因此改为手动求导
dg = ones(2*n, 1);
dg(n+1:2*n) = -1;
gxk(gxk>0) = 0;
dgxk = gxk'*dg;
y = f(x) - sum(lambda_h'*h(x)) + 0.5*mu*sum(h(x).^2);   % y=L(x);
grad=matlabFunction(diff(y));   %通过matlabFunction将符号函数转换为匿名函数，并求导
epsilonk = [grad(xk)+dgxk; h(xk); min(g(xk),lambda_g)];
epsilonk = norm(epsilonk);
t2=clock;
time = etime(t2,t1);
fprintf('\t问题编号\t\t\t算法\t\t运行次数\t\t运行时间\t\t\tKKT残量\n');
fprintf('\t%s\t\t乘子法\t%5d\t\t%2f\t\t%5d\n', nprob, iter_num, time, epsilonk);
return;
end

二、 Frank-Wolfe 算法

简介
解约束问题可行方向法与求解无约束问题的下降算法类似，但对约束问题而言，我们关心的是问题的可行点。可行方向法从问题的可行点出发，在该点的可行方向中，寻找使目标函数下降的方向。然后沿该方向进行线性搜索，得到一个可行点。如此得到一个点列，并希望该点列能终止于问题的解，或其极限点是问题的解。
Frank-Wolfe 算法便是利用求解约束问题可行方向法的思想，该算法通过解下面线性规划问题来确定可行下降方向。其搜索方向d(k)d^{(k)}d(k)可理解为可行点x(k)x^{(k)}x(k)处的可行方向中的最速下降法，是最速下降法的一种推广。
min⁡y∈DLf(x(k))+∇f(x(k))T(y−x(k))\min_{y \in D_L} f(x^{(k)}) + \nabla f(x^{(k)})^T(y - x^{(k)})y∈DLminf(x(k))+∇f(x(k))T(y−x(k))
或等价的
min⁡y∈DL∇f(x(k))T(y−x(k))\min_{y \in D_L} \nabla f(x^{(k)})^T(y - x^{(k)})y∈DLmin∇f(x(k))T(y−x(k))
其中DLD_LDL是问题的可行域.
适用范围：适用于求解约束函数为线性函数的最优化问题
算法步骤
步0：取初始点x(0)∈DLx^{(0)} \in D_Lx(0)∈DL, 精度ϵ>0\epsilon > 0ϵ>0. 令k=0.
步1：求解下面线性规划问题得y(k)y^{(k)}y(k). 令d(k)=y(k)−x(k)d^{(k)} = y^{(k)} - x^{(k)}d(k)=y(k)−x(k).
min⁡y∈DL∇f(x(k))T(y−x(k))\min_{y \in D_L} \nabla f(x^{(k)})^T(y - x^{(k)})y∈DLmin∇f(x(k))T(y−x(k))
步2：若∣∇f(x(k))Td(k)∣≤ϵ|\nabla f(x^{(k)})^Td^{(k)}| \leq \epsilon∣∇f(x(k))Td(k)∣≤ϵ，则得解x(k)x^{(k)}x(k). 否则，转步3.
步3：由线性搜索
min⁡0≤t≤1f(x(k)+td(k))=f(x(k)+tkd(k))\min_{0\leq t \leq 1} f(x^{(k)} + td^{(k)}) = f(x^{(k)} + t_kd^{(k)})0≤t≤1minf(x(k)+td(k))=f(x(k)+tkd(k))确定步长tkt_ktk. 令x(k+1)=x(k)+tkd(k),k=k+1x^{(k+1)} = x^{(k)} + t_kd^{(k)}, k=k+1x(k+1)=x(k)+tkd(k),k=k+1. 转步1.
matlab实现
说明：函数中使用的initpt函数（用于初始化x0x_0x0）是为了自己编写的辅助函数，这几个函数的具体实现可见最下方的链接处。

function [iter_num, time, epsilonk] = Frank_Wolfe(nprob, n, x, epsilon, max_iter)
% Frank Wolfe算法   2020.06.25
% @author: 豆奶
% 函数功能：实现Frank Wolfe法
% 输入：
% nprob: 问题编号，
% n: 维度
% x: 自变量x
% epsilon: int,
% max_iter: 最大迭代次数
% 输出：
% iter_num: 迭代次数
% time: 计算时间
% epsilonk: KKT系统模的最终值t1=clock;   % 计时开始
if nargin==2% 步0：取初始点x和精度epsilonx = initpt(n, nprob);  % 选取初始点要注意在可行域内 epsilon = 0.0001;max_iter = 100;
end
xk = x;
% 设置目标函数，约束函数，这一部分f和df也可以改成参数的形式传入.
if strcmp(nprob, 'De_Jong')f = @(x) sum(x.^2);df = @(x) 2.*x;g1 =@(x) x+5.12;g2 =@(x) 5.12-x;g =@(x) [g1(x); g2(x)];h =@(x) 0;
elseif strcmp(nprob, 'Rastrigin')f = @(x)10*n + sum(x.^2 - 10*cos(2*pi*x));   % n为常量df = @(x) 2*x + 20*pi*sin(2*pi*x);   % f 的导数g1 =@(x) 5.12+x;g2 =@(x) 5.12-x;g =@(x) [g1(x); g2(x)];h =@(x) 0;
end
syms t;
for k=0:max_iterdfx = df(xk);lb = -5.12*ones(n, 1);ub = 5.12*ones(n, 1);% 步1：利用matlab的linprog函数求解线性规划问题（13.8）得y% 参数说明： linprog(f,A,b,Aeq,beq,LB,UB, X0)  y = linprog(dfx,[],[],[],[],lb,ub);% 这里要改（6.22）, 改好了（6.22）d = y - xk;% 步2：判断是否满足终止条件if norm(dfx'*d)/(1+norm(f(xk)))<epsilonbreak;end% 步3：利用matlab的fmincon函数求解线性搜索确定步长tft = matlabFunction(f(xk+t*d));  %通过matlabFunction将符号函数转换为匿名函数，因为fmincon的第一个输入是匿名函数% 参数说明： fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options) %最后两个输入这里不需要[t_ans] = fmincon(ft,0.5,[],[],[],[],0,1);  % 这里要改（6.21）  改好了（6.22）xk = xk + t_ans*d;
end
% 获取输出值
iter_num = k;
lambda_g = zeros(2*n, 1);  % 由于Frank_Wolfe算法没有用到lagrange函数，故lambda设为0
epsilonk = [df(xk); h(xk); min(g(xk), lambda_g)];
epsilonk = norm(epsilonk);
t2=clock;
time = etime(t2,t1);
fprintf('\t问题名称\t\t\t算法\t\t运行次数\t\t运行时间\t\t\tKKT系统残量\n');
fprintf('\t%s\t\tFrank_Wolfe\t%5d\t\t%2f\t\t%5d\n', nprob, iter_num, time, epsilonk);
return;
end

参考文献：

数值最优化算法与理论(第二版)

完整代码与函数问题测试github(点我)

如有不对之处，还望指出，我会进行修改的，谢谢！

笔记：常见的约束问题求解算法——乘子法和Frank-Wolfe算法相关推荐

【数据结构笔记35】C实现：有序子列的归并算法：递归与非递归的实现
本次笔记内容: 9.4.1 有序子列的归并 9.4.2 归并算法 9.4.3 非递归算法文章目录归并排序的核心:两个有序子列的归并归并算法实现策略一:递归算法注意要有统一函数接口为什么要有参 ...
01背包问题的动态规划算法、蛮力法和空间优化算法
算法思想: (1).动态规划算法:解决背包物品价值最大化问题的最优解,是建立在每一个子问题的最优解的前提下完成的.设Value[i,j]表示的是i个物品放进背包容量为j的背包的价值,令i从0增至n(物 ...
直观理解拉格朗日乘子法和Karush-Kuhn-Tucker(KKT)条件
在最优化问题中,经常是会有约束条件的,而约束条件可分为等式约束条件和不等式约束条件,对于前者,我们有拉格朗日乘子法,对于后者,有KKT条件,对于既有等式约束又有不等式约束的最优化问题,只需要结合拉格朗 ...
学习笔记--一个自管理（组织）的多目标进化算法（SMEA）
学习笔记–一个自管理(组织)的多目标进化算法(SMEA) 摘要:在温和条件下,一个连续m维目标的优化问题的帕累托前沿(解集)可以形成一个(m-1)维的分段连续流形.基于这个性质,这篇文章提出了一个自管 ...
MySQL学习笔记06【多表查询、子查询、多表查询练习】
MySQL 文档-黑马程序员(腾讯微云):https://share.weiyun.com/RaCdIwas 1-MySQL基础.pdf.2-MySQL约束与设计.pdf.3-MySQL多表查询与事务 ...
算法基础：图的相关算法知识笔记
一.图的相关算法 1.图的分类知识如下图: 2.生成树概念对连通图进行遍历,过程中所经过的边和顶点的组合可看做是一棵普通树,通常称为生成树. 连通图的生成树具有这样的特征:边的数量 = 顶点数 - ...
算法基础：常用的排序算法知识笔记
1.算法外排序分类 2.冒泡排序冒泡排序(Bubble Sort)属于交换排序,它的原理是:循环两两比较相邻的记录,如果反序则交换,直到没有反序的记录为止. 实现算法: /** * 冒泡排序优化后的 ...
算法查找子节点_掌握着十大编程算法助你走上高手之路
算法一:快速排序算法快速排序是由东尼·霍尔所发展的一种排序算法.在平均状况下,排序 n 个项目要 Ο(n log n)次比较.在最坏状况下则需要 Ο(n2)次比较,但这种状况并不常见.事实上,快速排 ...
从零开始学习VIO笔记 --- 第四讲：滑动窗口（基于滑动窗口算法的 VIO 系统：可观性和一致性）
从零开始学习VIO笔记 --- 第四讲:滑动窗口(基于滑动窗口算法的 VIO 系统:可观性和一致性) 一. 从高斯分布到信息矩阵 1.1 高斯分布 1.2 高斯分布和协方差矩阵 1.3 信息矩阵二. ...
分类算法列一下有多少种？应用场景？分类算法介绍、常见分类算法优缺点、如何选择分类算法、分类算法评估
分类算法分类算法介绍概念分类算法常见分类算法 NBS LR SVM算法 ID3算法 C4.5 算法 C5.0算法 KNN 算法 ANN 算法选择分类算法分类算法性能评估分类算法介绍概念 ...

笔记：常见的约束问题求解算法——乘子法和Frank-Wolfe算法

笔记：常见的约束问题求解算法——乘子法和Frank-Wolfe算法相关推荐

最新文章

热门文章