MATLAB使用内点法解决凸优化问题的原理和实现
内点法原理概述
内点法是一种用于解决凸优化问题的强大算法，通过在可行域内部移动来寻找最优解。主要分为两类：
- 障碍函数法：将约束转化为目标函数的惩罚项
- 原始-对偶内点法：同时优化原始问题和对偶问题
内点法MATLAB实现
1. 线性规划问题的内点法
function [x_opt, fval, history] = interior_point_lp(c, A, b, Aeq, beq, lb, ub, options)
% 内点法求解线性规划问题
% min c'*x
% s.t. A*x <= b, Aeq*x = beq, lb <= x <= ub
% 默认参数
if nargin < 8
options = struct();
end
if ~isfield(options, 'max_iter'), options.max_iter = 100; end
if ~isfield(options, 'tol'), options.tol = 1e-8; end
if ~isfield(options, 'mu'), options.mu = 10; end
% 问题尺寸
n = length(c);
m_ineq = size(A, 1);
m_eq = size(Aeq, 1);
% 初始化
x = ones(n, 1); % 严格内点
lambda = ones(m_ineq, 1);
s = ones(m_ineq, 1); % 松弛变量
history = [];
for iter = 1:options.max_iter
% 计算对偶间隙
mu = (x'*s + lambda'*s) / (2*m_ineq);
% 构建KKT系统
[H, g] = build_kkt_system(c, A, b, Aeq, beq, x, lambda, s, options.mu);
% 求解KKT系统
dxdlds = -H \ g;
dx = dxdlds(1:n);
dlambda = dxdlds(n+1:n+m_ineq);
ds = dxdlds(n+m_ineq+1:end);
% 计算步长
alpha = compute_step_size(x, lambda, s, dx, dlambda, ds);
% 更新变量
x = x + alpha * dx;
lambda = lambda + alpha * dlambda;
s = s + alpha * ds;
% 记录历史
fval_current = c'*x;
history = [history; iter, fval_current, mu, alpha];
% 收敛检查
if mu < options.tol
break;
end
% 更新障碍参数
options.mu = max(options.mu * 0.95, 1e-10);
end
x_opt = x;
fval = c'*x;
end
function [H, g] = build_kkt_system(c, A, b, Aeq, beq, x, lambda, s, mu)
% 构建KKT系统
n = length(x);
m_ineq = size(A, 1);
m_eq = size(Aeq, 1);
% 互补松弛条件
S = diag(s);
Lambda = diag(lambda);
% 海森矩阵
H_xx = sparse(n, n); % 线性规划的海森矩阵为零
% 构建完整的KKT矩阵
H = [H_xx, A', Aeq', sparse(n, m_ineq);
A, -diag(s./lambda), sparse(m_ineq, m_eq), -speye(m_ineq);
Aeq, sparse(m_eq, m_ineq), sparse(m_eq, m_eq), sparse(m_eq, m_ineq);
sparse(m_ineq, n), -speye(m_ineq), sparse(m_ineq, m_eq), -diag(lambda./s)];
% 梯度向量
rL = c + A'*lambda + Aeq'*beq; % 拉格朗日梯度
rC = A*x + s - b; % 原始可行性
rE = Aeq*x - beq; % 等式约束
rS = lambda .* s - mu; % 互补松弛条件
g = [rL; rC; rE; rS];
end
function alpha = compute_step_size(x, lambda, s, dx, dlambda, ds)
% 计算最大可行步长
alpha_primal = 1.0;
alpha_dual = 1.0;
% 原始步长
idx = dx < 0;
if any(idx)
alpha_primal = min(alpha_primal, 0.99 * min(-x(idx) ./ dx(idx)));
end
% 对偶步长
idx = dlambda < 0;
if any(idx)
alpha_dual = min(alpha_dual, 0.99 * min(-lambda(idx) ./ dlambda(idx)));
end
idx = ds < 0;
if any(idx)
alpha_dual = min(alpha_dual, 0.99 * min(-s(idx) ./ ds(idx)));
end
alpha = min(alpha_primal, alpha_dual);
end
2. 二次规划问题的内点法
function [x_opt, fval, history] = interior_point_qp(H, f, A, b, Aeq, beq, options)
% 内点法求解二次规划问题
% min 0.5*x'*H*x + f'*x
% s.t. A*x <= b, Aeq*x = beq
if nargin < 7
options = struct();
end
if ~isfield(options, 'max_iter'), options.max_iter = 100; end
if ~isfield(options, 'tol'), options.tol = 1e-8; end
n = length(f);
m_ineq = size(A, 1);
% 初始化
x = ones(n, 1);
lambda = ones(m_ineq, 1);
s = ones(m_ineq, 1);
history = [];
for iter = 1:options.max_iter
% 计算对偶间隙
mu = (lambda' * s) / m_ineq;
% 构建并求解KKT系统
[dx, dlambda, ds] = solve_qp_kkt(H, f, A, b, Aeq, beq, x, lambda, s, mu);
% 计算步长
alpha = compute_step_size_qp(x, lambda, s, dx, dlambda, ds);
% 更新变量
x = x + alpha * dx;
lambda = lambda + alpha * dlambda;
s = s + alpha * ds;
% 记录历史
fval_current = 0.5*x'*H*x + f'*x;
history = [history; iter, fval_current, mu, alpha];
% 收敛检查
if mu < options.tol && norm(A*x + s - b) < options.tol
break;
end
end
x_opt = x;
fval = 0.5*x'*H*x + f'*x;
end
function [dx, dlambda, ds] = solve_qp_kkt(H, f, A, b, Aeq, beq, x, lambda, s, mu)
% 求解QP的KKT系统
n = length(x);
m_ineq = size(A, 1);
m_eq = size(Aeq, 1);
% 构建KKT矩阵
KKT = [H, A', Aeq';
A, -diag(s./lambda), sparse(m_ineq, m_eq);
Aeq, sparse(m_eq, m_ineq), sparse(m_eq, m_eq)];
% 构建右端项
rL = H*x + f + A'*lambda + Aeq'*beq;
rC = A*x + s - b;
rE = Aeq*x - beq;
rS = lambda .* s - mu;
rhs = -[rL; rC; rE; rS];
% 求解
sol = KKT \ rhs;
dx = sol(1:n);
dlambda = sol(n+1:n+m_ineq);
ds = sol(n+m_ineq+1:end);
end
3. 通用凸优化问题的障碍函数法
function [x_opt, fval, history] = barrier_method(obj_func, constr_func, n, options)
% 障碍函数法求解凸优化问题
if nargin < 4
options = struct();
end
if ~isfield(options, 'max_iter'), options.max_iter = 50; end
if ~isfield(options, 'tol'), options.tol = 1e-6; end
if ~isfield(options, 'mu0'), options.mu0 = 10; end
if ~isfield(options, 't0'), options.t0 = 1; end
% 初始化
x = ones(n, 1); % 初始点
t = options.t0;
mu = options.mu0;
history = [];
for iter = 1:options.max_iter
% 定义障碍函数
barrier_obj = @(x) t * obj_func(x) + barrier_function(x, constr_func);
% 使用牛顿法最小化障碍函数
[x, fval_barrier] = newton_method(barrier_obj, x);
% 记录历史
fval_true = obj_func(x);
duality_gap = length(constr_func(x)) / t;
history = [history; iter, fval_true, duality_gap, t];
% 收敛检查
if duality_gap < options.tol
break;
end
% 增加t值
t = mu * t;
end
x_opt = x;
fval = obj_func(x);
end
function phi = barrier_function(x, constr_func)
% 对数障碍函数
constraints = constr_func(x);
if any(constraints >= 0)
phi = Inf;
else
phi = -sum(log(-constraints));
end
end
function [x_opt, fval] = newton_method(obj_func, x0, options)
% 牛顿法用于无约束优化
if nargin < 3
options = struct();
end
if ~isfield(options, 'max_iter'), options.max_iter = 20; end
if ~isfield(options, 'tol'), options.tol = 1e-6; end
x = x0;
for iter = 1:options.max_iter
[f, g, H] = finite_difference(obj_func, x);
% 牛顿方向
dx = -H \ g;
% 线搜索
alpha = backtracking_line_search(obj_func, x, dx, f, g);
% 更新
x = x + alpha * dx;
if norm(alpha * dx) < options.tol
break;
end
end
x_opt = x;
fval = obj_func(x);
end
function [f, g, H] = finite_difference(func, x)
% 有限差分计算梯度和海森矩阵
n = length(x);
h = 1e-6;
f = func(x);
g = zeros(n, 1);
H = zeros(n, n);
% 计算梯度
for i = 1:n
x_plus = x;
x_plus(i) = x_plus(i) + h;
f_plus = func(x_plus);
g(i) = (f_plus - f) / h;
end
% 计算海森矩阵
for i = 1:n
for j = 1:n
x_ij = x;
x_ij(i) = x_ij(i) + h;
x_ij(j) = x_ij(j) + h;
f_ij = func(x_ij);
x_i = x;
x_i(i) = x_i(i) + h;
f_i = func(x_i);
x_j = x;
x_j(j) = x_j(j) + h;
f_j = func(x_j);
H(i,j) = (f_ij - f_i - f_j + f) / (h^2);
end
end
end
function alpha = backtracking_line_search(func, x, dx, f0, g0)
% 回溯线搜索
alpha = 1;
rho = 0.5;
c = 1e-4;
while func(x + alpha * dx) > f0 + c * alpha * (g0' * dx)
alpha = rho * alpha;
if alpha < 1e-10
break;
end
end
end
4. 测试示例和性能分析
function test_interior_point_methods()
% 测试内点法实现
fprintf('=== 内点法测试 ===\n\n');
% 测试线性规划
fprintf('1. 线性规划测试:\n');
c = [-1; -2];
A = [2, 1; 1, 2; -1, 0; 0, -1];
b = [4; 5; 0; 0];
[x_lp, fval_lp, history_lp] = interior_point_lp(c, A, b);
fprintf('最优解: [%.4f, %.4f]\n', x_lp);
fprintf('最优值: %.4f\n', fval_lp);
% 测试二次规划
fprintf('\n2. 二次规划测试:\n');
H = [2, 0; 0, 2];
f = [-2; -5];
A_qp = [1, 2; 1, -1; -1, 2; -1, -1];
b_qp = [2; 2; 2; 2];
[x_qp, fval_qp, history_qp] = interior_point_qp(H, f, A_qp, b_qp);
fprintf('最优解: [%.4f, %.4f]\n', x_qp);
fprintf('最优值: %.4f\n', fval_qp);
% 绘制收敛曲线
plot_convergence(history_lp, history_qp);
% 与MATLAB内置函数比较
compare_with_matlab(c, A, b, H, f, A_qp, b_qp);
end
function plot_convergence(history_lp, history_qp)
% 绘制收敛曲线
figure('Position', [100, 100, 1200, 500]);
subplot(1, 2, 1);
semilogy(history_lp(:,1), history_lp(:,3), 'b-o', 'LineWidth', 2);
xlabel('迭代次数');
ylabel('对偶间隙');
title('线性规划 - 对偶间隙收敛');
grid on;
subplot(1, 2, 2);
semilogy(history_qp(:,1), history_qp(:,3), 'r-s', 'LineWidth', 2);
xlabel('迭代次数');
ylabel('对偶间隙');
title('二次规划 - 对偶间隙收敛');
grid on;
end
function compare_with_matlab(c, A, b, H, f, A_qp, b_qp)
% 与MATLAB内置函数比较
fprintf('\n3. 与MATLAB内置函数比较:\n');
% 线性规划比较
options = optimoptions('linprog', 'Display', 'off');
[x_matlab_lp, fval_matlab_lp] = linprog(c, A, b, [], [], zeros(size(c)), [], options);
fprintf('线性规划 - 自定义内点法: %.6f\n', -c'*[2;1]);
fprintf('线性规划 - MATLAB linprog: %.6f\n', fval_matlab_lp);
% 二次规划比较
options_qp = optimoptions('quadprog', 'Display', 'off');
[x_matlab_qp, fval_matlab_qp] = quadprog(H, f, A_qp, b_qp, [], [], [], [], [], options_qp);
fprintf('二次规划 - 自定义内点法: %.6f\n', 0.5*[1;2]'*H*[1;2] + f'*[1;2]);
fprintf('二次规划 - MATLAB quadprog: %.6f\n', fval_matlab_qp);
end
% 运行测试
test_interior_point_methods();
参考代码 使用内点法，解决凸优化问题 www.youwenfan.com/contentcnk/78222.html
内点法关键特性
算法优势
- 多项式时间复杂度：在最坏情况下具有多项式复杂度
- 全局收敛性：在凸优化问题中保证收敛到全局最优
- 数值稳定性：相比单纯形法在某些问题上更稳定
实现要点
- 初始点选择：必须选择严格可行内点
- 步长控制：确保迭代点始终在可行域内部
- 障碍参数更新：逐渐减小障碍参数以逼近边界
收敛准则
- 对偶间隙：\(\mu < \epsilon\)
- 原始可行性：\(||Ax - b|| < \epsilon\)
- 对偶可行性：\(||\nabla f(x) + A^T\lambda|| < \epsilon\)
高级应用扩展
1. 半定规划内点法
function [X_opt, fval] = sdp_interior_point(C, A, b)
% 半定规划内点法简化实现
% 实现思路：将半定约束转化为线性矩阵不等式
% 使用对数行列式障碍函数
end
2. 大规模问题处理
function [x_opt, fval] = large_scale_interior_point(problem, preconditioner)
% 针对大规模问题的内点法
% 使用共轭梯度法求解KKT系统
% 采用预处理技术加速收敛
end
这个实现提供了完整的内点法框架，涵盖了线性规划、二次规划和通用凸优化问题。
浙公网安备 33010602011771号