优化准则法和数学规划法
优化准则法和数学规划法，包括它们的原理、对比和应用示例
优化方法总体分类
首先通过一个表格来概览这两种方法的对比：
| 特征维度 | 优化准则法 | 数学规划法 |
|---|---|---|
| 基本思想 | 基于工程直觉和物理准则 | 基于数学理论和算法 |
| 收敛性 | 通常较快，但可能收敛到局部最优 | 收敛性有理论保证，可找到全局最优 |
| 适用问题 | 结构优化、拓扑优化等 | 各类优化问题，通用性强 |
| 实现难度 | 相对简单，物理意义明确 | 需要较强的数学基础 |
| 典型方法 | 满应力准则、能量准则等 | 梯度下降、序列二次规划等 |
优化准则法详解
1. 基本原理
优化准则法基于工程直觉和物理最优性条件，通过构造最优性准则来指导优化过程。
2. 满应力准则法示例
function [x_opt, history] = optimality_criterion_method()
% 满应力准则法 - 桁架结构优化
% 最小化重量，满足应力约束
% 问题参数
n_elements = 10; % 单元数量
max_iter = 100; % 最大迭代次数
tol = 1e-6; % 收敛容差
rho = 7800; % 材料密度 kg/m^3
L = 1; % 单元长度 m
sigma_max = 250e6; % 最大允许应力 Pa
sigma_min = -200e6; % 最小允许应力 Pa
% 初始化设计变量（截面面积）
x = ones(n_elements, 1) * 1e-3; % 初始截面面积 m^2
x_min = 1e-4 * ones(n_elements, 1); % 最小截面面积
x_max = 1e-2 * ones(n_elements, 1); % 最大截面面积
% 存储优化历史
history.x = zeros(n_elements, max_iter);
history.obj = zeros(1, max_iter);
history.max_stress = zeros(1, max_iter);
% 优化循环
for iter = 1:max_iter
% 1. 结构分析 - 计算单元应力
stress = structural_analysis(x);
% 2. 计算目标函数（总重量）
weight = rho * L * sum(x);
% 3. 应用满应力准则更新设计变量
x_new = zeros(n_elements, 1);
for i = 1:n_elements
if stress(i) > 0 % 拉应力
x_new(i) = x(i) * abs(stress(i)) / sigma_max;
else % 压应力
x_new(i) = x(i) * abs(stress(i)) / abs(sigma_min);
end
end
% 4. 应用尺寸约束
x_new = max(x_min, min(x_max, x_new));
% 5. 检查收敛
if norm(x_new - x) < tol
fprintf('收敛于第 %d 次迭代\n', iter);
break;
end
% 6. 阻尼更新（避免振荡）
alpha = 0.7; % 阻尼系数
x = alpha * x_new + (1 - alpha) * x;
% 存储历史
history.x(:, iter) = x;
history.obj(iter) = weight;
history.max_stress(iter) = max(abs(stress));
fprintf('迭代 %d: 重量 = %.4f kg, 最大应力 = %.2f MPa\n', ...
iter, weight, max(abs(stress))/1e6);
end
x_opt = x;
history.iter = iter;
end
function stress = structural_analysis(x)
% 简化的桁架结构分析
% 在实际应用中应使用有限元分析
n = length(x);
stress = zeros(n, 1);
% 模拟应力分布（简化的应力计算）
base_load = 1e5; % 基础载荷 N
for i = 1:n
% 应力 = 载荷 / 面积，加上一些变化
variation = 0.5 * sin(2*pi*i/n);
nominal_stress = base_load / x(i);
stress(i) = nominal_stress * (1 + variation);
end
end
3. 能量准则法示例
function compliance_minimization()
% 基于柔度最小化的优化准则法
% 问题设置
n_elements = 8;
max_iter = 50;
vol_frac = 0.4; % 体积分数约束
% 设计变量初始化（单元密度）
x = ones(n_elements, 1) * vol_frac;
x_min = 0.001;
% 优化参数
move = 0.2; % 移动限制
p = 3; % SIMP惩罚因子
for iter = 1:max_iter
% 计算柔度和灵敏度
[C, dC] = compute_compliance(x, p);
% 优化准则更新（OC方法）
l1 = 0; l2 = 1e9;
while (l2 - l1) / (l1 + l2) > 1e-4
lmid = 0.5 * (l1 + l2);
x_new = max(x_min, max(x - move, min(1, min(x + move, ...
x .* sqrt(-dC / lmid))));
if sum(x_new) > vol_frac * n_elements
l1 = lmid;
else
l2 = lmid;
end
end
% 更新设计变量
x = x_new;
fprintf('迭代 %d: 柔度 = %.6f, 体积分数 = %.3f\n', ...
iter, C, mean(x));
end
end
function [C, dC] = compute_compliance(x, p)
% 计算柔度和灵敏度（简化版本）
n = length(x);
% SIMP材料插值
E = x.^p; % 弹性模量
% 简化的柔度计算
C = sum(1 ./ E);
% 灵敏度分析
dC = -p * x.^(p-1) ./ E.^2;
end
数学规划法详解
1. 梯度下降法示例
function [x_opt, history] = gradient_descent()
% 梯度下降法求解无约束优化问题
% 最小化 Rosenbrock 函数
% 目标函数: f(x) = (1-x1)^2 + 100*(x2-x1^2)^2
f = @(x) (1-x(1))^2 + 100*(x(2)-x(1)^2)^2;
grad = @(x) [2*(x(1)-1) - 400*x(1)*(x(2)-x(1)^2);
200*(x(2)-x(1)^2)];
% 优化参数
max_iter = 1000;
tol = 1e-8;
x0 = [-1; 1]; % 初始点
% 线搜索参数
alpha0 = 1; % 初始步长
c = 1e-4; % Armijo条件参数
rho = 0.5; % 步长缩减因子
x = x0;
history.x = zeros(2, max_iter);
history.f = zeros(1, max_iter);
history.grad_norm = zeros(1, max_iter);
for iter = 1:max_iter
% 计算梯度和函数值
g = grad(x);
f_val = f(x);
% 存储历史
history.x(:, iter) = x;
history.f(iter) = f_val;
history.grad_norm(iter) = norm(g);
% 检查收敛
if norm(g) < tol
fprintf('梯度收敛于第 %d 次迭代\n', iter);
break;
end
% 线搜索确定步长
alpha = alpha0;
while f(x - alpha*g) > f_val - c*alpha*(g'*g)
alpha = rho * alpha;
end
% 更新变量
x = x - alpha * g;
if mod(iter, 100) == 0
fprintf('迭代 %d: f(x) = %.6f, ||∇f|| = %.6f\n', ...
iter, f_val, norm(g));
end
end
x_opt = x;
history.iter = iter;
end
2. 序列二次规划(SQP)示例
function [x_opt, history] = sqp_method()
% 序列二次规划法求解约束优化问题
% 问题定义
% 最小化 f(x) = x1^2 + x2^2
% 约束: x1 + x2 >= 1
% x1^2 + x2^2 <= 4
f = @(x) x(1)^2 + x(2)^2;
grad_f = @(x) [2*x(1); 2*x(2)];
hess_f = @(x) [2, 0; 0, 2];
% 约束函数
g1 = @(x) x(1) + x(2) - 1; % g1(x) >= 0
g2 = @(x) 4 - x(1)^2 - x(2)^2; % g2(x) >= 0
grad_g1 = @(x) [1; 1];
grad_g2 = @(x) [-2*x(1); -2*x(2)];
% SQP参数
max_iter = 50;
tol = 1e-6;
x = [0; 0]; % 初始点
lambda = [0; 0]; % 拉格朗日乘子
history.x = zeros(2, max_iter);
history.f = zeros(1, max_iter);
history.violation = zeros(1, max_iter);
for iter = 1:max_iter
% 当前点的梯度和Hessian
g_f = grad_f(x);
H = hess_f(x);
% 约束函数值和梯度
g_val = [g1(x); g2(x)];
A = [grad_g1(x)'; grad_g2(x)'];
% 求解QP子问题
% min 0.5*d'*H*d + g_f'*d
% s.t. A*d >= -g_val
options = optimoptions('quadprog', 'Display', 'none');
d = quadprog(H, g_f, -A, g_val, [], [], [], [], [], options);
% 更新变量
x_new = x + d;
% 更新拉格朗日乘子
lambda_new = -A \ g_f;
% 检查收敛
if norm(x_new - x) < tol && norm(lambda_new - lambda) < tol
fprintf('SQP收敛于第 %d 次迭代\n', iter);
break;
end
x = x_new;
lambda = lambda_new;
% 存储历史
history.x(:, iter) = x;
history.f(iter) = f(x);
history.violation(iter) = max(-min(g_val), 0);
fprintf('迭代 %d: f(x) = %.4f, 约束违反 = %.6f\n', ...
iter, f(x), history.violation(iter));
end
x_opt = x;
history.iter = iter;
end
3. 内点法示例
function [x_opt, history] = interior_point_method()
% 内点法求解不等式约束优化问题
% 问题: 最小化 f(x) = x1^4 + x2^4
% 约束: x1 + x2 >= 1
% x1 >= 0, x2 >= 0
f = @(x) x(1)^4 + x(2)^4;
grad_f = @(x) [4*x(1)^3; 4*x(2)^3];
hess_f = @(x) [12*x(1)^2, 0; 0, 12*x(2)^2];
% 内点法参数
mu0 = 10; % 初始障碍参数
mu = mu0;
sigma = 0.1; % 缩减因子
tol = 1e-6;
max_iter = 100;
x = [0.5; 0.5]; % 严格内点
history.x = zeros(2, max_iter);
history.f = zeros(1, max_iter);
history.mu = zeros(1, max_iter);
for iter = 1:max_iter
% 障碍函数: f(x) - mu * sum(log(约束))
barrier_obj = @(x) f(x) - mu * (log(x(1)) + log(x(2)) + log(x(1)+x(2)-1));
barrier_grad = @(x) grad_f(x) - mu * [1/x(1) + 1/(x(1)+x(2)-1);
1/x(2) + 1/(x(1)+x(2)-1)];
barrier_hess = @(x) hess_f(x) + mu * [1/x(1)^2+1/(x(1)+x(2)-1)^2, 1/(x(1)+x(2)-1)^2;
1/(x(1)+x(2)-1)^2, 1/x(2)^2+1/(x(1)+x(2)-1)^2];
% 求解无约束子问题（使用牛顿法）
x_new = newton_method(barrier_obj, barrier_grad, barrier_hess, x);
% 更新障碍参数
mu = sigma * mu;
% 检查收敛
if norm(x_new - x) < tol
fprintf('内点法收敛于第 %d 次迭代\n', iter);
break;
end
x = x_new;
% 存储历史
history.x(:, iter) = x;
history.f(iter) = f(x);
history.mu(iter) = mu;
fprintf('迭代 %d: f(x) = %.6f, mu = %.6f\n', iter, f(x), mu);
end
x_opt = x;
history.iter = iter;
end
function x_opt = newton_method(obj, grad, hess, x0)
% 牛顿法求解无约束优化问题
max_iter = 20;
tol = 1e-8;
x = x0;
for i = 1:max_iter
g = grad(x);
H = hess(x);
if norm(g) < tol
break;
end
d = -H \ g; % 牛顿方向
% 线搜索
alpha = 1;
while obj(x + alpha*d) > obj(x) + 1e-4*alpha*g'*d
alpha = 0.5 * alpha;
end
x = x + alpha * d;
end
x_opt = x;
end
综合对比与可视化
function compare_methods()
% 比较不同优化方法的性能
% 测试函数: f(x) = x1^2 + x2^2 + x1*x2
f = @(x) x(1)^2 + x(2)^2 + x(1)*x(2);
grad_f = @(x) [2*x(1) + x(2); 2*x(2) + x(1)];
% 不同初始点
start_points = [-2, 2; 2, -2; -1, -1];
figure('Position', [100, 100, 1200, 800]);
% 绘制等高线
subplot(2, 3, 1);
[X, Y] = meshgrid(-3:0.1:3, -3:0.1:3);
Z = arrayfun(@(x,y) f([x,y]), X, Y);
contour(X, Y, Z, 50);
hold on;
title('目标函数等高线');
xlabel('x1'); ylabel('x2');
% 测试不同方法
methods = {'梯度下降', '牛顿法', '共轭梯度'};
for i = 1:size(start_points, 1)
x0 = start_points(i, :)';
for method_idx = 1:3
switch method_idx
case 1
[x_opt, hist] = run_gradient_descent(f, grad_f, x0);
color = 'r';
case 2
[x_opt, hist] = run_newton_method(f, grad_f, x0);
color = 'g';
case 3
[x_opt, hist] = run_conjugate_gradient(f, grad_f, x0);
color = 'b';
end
% 绘制优化路径
subplot(2, 3, 1);
plot(hist.x(1, 1:hist.iter), hist.x(2, 1:hist.iter), ...
[color, '-o'], 'LineWidth', 1, 'MarkerSize', 3);
% 绘制收敛曲线
subplot(2, 3, method_idx + 1);
semilogy(1:hist.iter, hist.obj(1:hist.iter), ...
[color(1), '-'], 'LineWidth', 2);
hold on;
title([methods{method_idx}, ' 收敛曲线']);
xlabel('迭代次数'); ylabel('目标函数值');
grid on;
end
end
subplot(2, 3, 1);
legend('等高线', 'GD-路径1', 'Newton-路径1', 'CG-路径1', ...
'GD-路径2', 'Newton-路径2', 'CG-路径2', ...
'GD-路径3', 'Newton-路径3', 'CG-路径3');
end
参考代码 优化准则法和数学规划法 www.youwenfan.com/contentcnm/79812.html
应用建议
-
选择优化准则法的情况：
- 问题具有明确的物理意义和最优性准则
- 需要快速获得可行解
- 工程经验丰富的领域（如结构优化）
-
选择数学规划法的情况：
- 问题具有复杂的数学约束
- 需要理论上的收敛保证
- 通用性要求较高的问题
-
混合方法：
- 使用优化准则法获得好的初始点
- 用数学规划法进行精细优化
- 结合两种方法的优点
浙公网安备 33010602011771号