多目标规划问题求解
全糖奶茶屋
2022年09月07日 18:16

        多目标规划问题, 就是在规划问题的基础上, 增加为多个目标函数, 比如除了市场经营使得盈利最大的同时, 保证风险最小. 同规划模型一样, 它的建模过程的关键, 依旧是决策变量, 目标函数和约束条件三部分, 写法也可以参照规划模型的写作手法来写, 最终在汇总之后, 我们会得到如下这样的一个标准形式:

多目标规划模型的标准形式

        接下来, 就需要介绍几种方法, 对上面的标准形式进行求解:

1. 传统解法——多目标规划转换为单目标规划

① 理想点法

        多目标规划的每个目标函数都有最优解，我们记录每个目标函数的最优解，每个目标函数的最优解构成的向量称为理想解，重新构造一个目标函数，这个目标函数的目的是使自变量的取值在可行解的范围内，距离理想解最近。

        以如下例题为例进行讲解:

例1: 利用理想点法求解规划模型:

规划模型例题

解: 先对单目标进行求解

        对第一个目标求解的matlab程序为：

代码块
JavaScript
自动换行
复制代码
f=[3;-2];
a=[2  3;2 1];
b=[18;10];
lb=[0;0];
ub=[];
[x,favl]=linprog(f,a,b,[],[],lb,ub);
复制成功

输出结果为：x=0.0000 6.0000，favl=-12.0000，所以最大值为12.0000

        对第二个目标求解的matlab程序为：

代码块
JavaScript
自动换行
复制代码
f=[-4;-3];
a=[2  3;2 1];
b=[18;10];
lb=[0;0];
ub=[];
[x,favl]=linprog(f,a,b,[],[],lb,ub);
复制成功

输出结果为：x=3.0000 4.0000，favl=-24.0000，所以最大值为24.0000

        于是得到理想点：（12,24）, 然后求如下模型的最优解:

理想点法简化多目标规划

该步的Matlab程序为:

代码块
JavaScript
自动换行
复制代码
x0=[1;1];
a=[2 3;2 1];
b=[18;10];
lb=[0;0];
ub=[];
x=fmincon('((-3*x(1)+2*x(2)-12)^2+(4*x(1)+3*x(2)-24)^2)^(1/2)',x0,a,b,[],[],lb,ub);
f1=-3*x(1)+2*x(2);
f2=4*x(1)+3*x(2);
复制成功

 ② 线性加权法：

        在具有多个指标的问题中，人们总希望对那些相对重要的指标给予较大的权系数，因而将多目标向量问题转化为所有目标的加权求和的标量问题，基于这个现实，构造如下评价函数，即

min_%7Bx%20%5Cin%20D%7D%20Z(x)%3D%5Csum_%7Bi%3D1%7D%5Er%20%5Comega_i%20Z_i(x)

将它的最优解 x 作为原模型在线性加权和意义下的“最优解”。（ %5Comega_i为加权因子，其选取的方法很多，有专家打分法、容限法和加权因子分解法等）.

        实际上, 如果学习过评价模型的话, 可以发现, 这里是借助了一种评价模型的思想, 将多个目标经过评价过程, 增加权重, 变成一个目标函数.

例2∶对例1进行线性加权和法求解。（权系数分别取 ）

解∶构造如下评价函数，即求如下模型的最优解。

线性加权和法简化多目标规划问题

该过程的Matlab程序为:

代码块
JavaScript
自动换行
复制代码
f=[-0.5;-2.5];
a=[2 3;2 1];
b=[18;10];
lb=[0;0];
ub=[];
x=linprog(f,a,b,[],[],lb,ub);
复制成功

输出结果为：x1=0.0000,x2=6.0000,对应的目标值为f1=12.0000,f2=18.0000.

 ③ 最大最小法：

        在决策的时候，采取保守策略是稳妥的，即在最坏的情况下，寻求最好的结果，按照此想法，可以构造如下评价函数，即

%5Cphi(Z)%3D%5Cmax_%7B1%20%5Cleq%20i%20%5Cleq%20r%7D%20Z_i

然后求解:

%5Cmin_%7Bx%20%5Cin%20D%7D%20%5Cphi%5BZ(x)%5D%3D%5Cmin_%7Bx%20%5Cin%20D%7D%5Cmax_%7B1%20%5Cleq%20i%20%5Cleq%20r%7D%20Z_i(x)

并将它的最优解 x 作为原模型在最大最小意义下的"最优解&#​34;。

        仍然以shangshu 例题为例进行讲解:

例3: 利用最大最小法求解规划模型:

规划模型例题

我们直接使用代码实现:

        首先编写M函数文件：

代码块
JavaScript
自动换行
复制代码
function f=mutiplesubjiect(x)
f(1)=3*x(1)-2*x(2);
f(2)=-4*x(1)-3*x(2);
复制成功

然后输入

代码块
JavaScript
自动换行
复制代码
>> x0=[0;0];
a=[2 3;2 1];
b=[18;10];
lb=[0;0];
[x,favl]=fminimax('mutiplesubjiect',x0,a,b,[],[],lb,[])
复制成功

④约束条件法

        理论依据 ：若规划问题的某一目标可以给出一个可供选择的范围，则该目标就可以作为约束条件而被排除出目标组，进入约束条件组中。         假如，除第一个目标外，其余目标都可以提出一个可供选择的范围，则该多目标规划问题就可以转化为单目标规划问题：

若原目标函数为:

%5Cmax%20Z%3D(f_1%2Cf_2%2C%5Ccdots%2Cf_n)

则可以保留 f_1 为目标函数, 其他目标函数转化为约束条件

f_j%5E%5Cmin%20%5Cleq%20f_j%20%5Cleq%20f_j%5E%5Cmax

并且不断下调他们的边界, 保证其中有解.

2. 智能优化算法求解

基于精英策略的快速非支配排序遗传算法（NSGA-II）

这里的理论部分可以参照遗传算法的理论:

遗传算法课程视频链接​

这里直接给出对应的代码分析与例题:

清除所有变量（非必须，但注意变量不能和下面所用的冲突）

代码块
JavaScript
自动换行
复制代码
clear all;
clc;
复制成功

需求解模型的参数设置部分：（模型导入）

代码块
JavaScript
自动换行
复制代码
%% 模型设置
% 适应度函数的函数句柄
fitnessfcn=@Fun;
% 变量个数
nvars=4;
% 约束条件形式1：下限与上限（若无取空数组[]）
% lb<= X <= ub
lb=[0,0,0,0];
ub=[];

% 约束条件形式2：线性不等式约束（若无取空数组[]）
% A*X <= b 
A = [0    0 1 1
    -1/3  0 0 0
     0 -1/2 0 0
     0    0 0 0];

b = [48 ; 30 ; 30 ; 0];

% 约束条件形式3：线性等式约束（若无取空数组[]）
% Aeq*X == beq
Aeq=[1 1 0 0;0 0 0 0; 0 0 0 0; 0 0 0 0];
beq=[120;0;0;0];
复制成功

目标函数: (这一段需放在脚本最后或单独放在一个文件里)

代码块
JavaScript
自动换行
复制代码
function y=Fun(x)
	% y是目标函数向量。有几个目标函数y就有多少个维度（数组y的长度）
	% 因为gamultiobj是以目标函数分量取极小值为目标，
	% 因此有些取极大值的目标函数注意取相反数
	y(1)=-(x(1)*100/3 + x(3)*90/3  + x(2)*80/2+x(4)*70/2);
	y(2)=x(3)+x(4);
end
复制成功

求解器设置部分：

代码块
JavaScript
自动换行
复制代码
%% 求解器设置
% 最优个体系数paretoFraction
% 种群大小populationsize
% 最大进化代数generations
% 停止代数stallGenLimit
% 适应度函数偏差TolFun
% 函数gaplotpareto：绘制Pareto前沿 
options=gaoptimset('paretoFraction',0.3,'populationsize',200,'generations',300,'stallGenLimit',200,'TolFun',1e-10,'PlotFcns',@gaplotpareto);
复制成功

求解与结果输出部分：

代码块
JavaScript
自动换行
复制代码
%% 主求解
[x,fval]=gamultiobj(fitnessfcn,nvars,A,b,Aeq,beq,lb,ub,options)

%% 结果提取
% 因为gamultiobj是以目标函数分量取极小值为目标，
% 因此在y=Fun(x)里取相反数的目标函数再取相反数画出原始情况
plot(-fval(:,1),fval(:,2),'pr')
xlabel('f_1(x)')
ylabel('f_2(x)')
title('Pareto front')
grid on
复制成功

        上面程序编写完成后进行运行:

求解时间受求解器设置影响，可能会较长，请耐心等待

求解过程中会实时显示当前种群的情况：

当前种群示意图

        如果已经达到满意，也可点击stop按钮提前结束求解

        最后的求解结果，即Pareto最优解集储存在[x,fval]中，fval是x对应的目标函数值。fval大致构成了一条空间曲线——Pareto前沿。若各个解较为均匀分布，说明该图包含了大部分最优解情况，全局性优，适用性强。在满足Pareto最优的条件下，是没有办法在不让某一优化目标受损的情况下，令另一方目标获得更优的。所以这些解均为最优，对最优解的具体选择可以根据实际情况。

例题:

完整代码:

代码块
JavaScript
自动换行
复制代码
clear
clc
fitnessfcn=@Fun;
% 变量个数
nvars=4;
% lb<= X <= ub
lb=[0,0,0,0];
ub=[270 240 460 130];
% A*X <= b 
A = [[13 13.5 14 11.5]
    -[13 13.5 14 11.5]
      0 0 -1 0
     [0.015 0.02 0.018 0.011]];

b = [300*48 ; -300*40 ; -150 ; 20];

% Aeq*X = beq
Aeq=[];beq=[];
%最优个体系数paretoFraction
%种群大小populationsize
%最大进化代数generations
%停止代数stallGenLimit
%适应度函数偏差TolFun
%函数gaplotpareto：绘制Pareto前沿 
options=gaoptimset('paretoFraction',0.3,'populationsize',200,'generations',300,'stallGenLimit',200,'TolFun',1e-10,'PlotFcns',@gaplotpareto);

[x,fval]=gamultiobj(fitnessfcn,nvars,A,b,Aeq,beq,lb,ub,options)

plot(-fval(:,1),fval(:,2),'pr')
xlabel('f_1(x)')
ylabel('f_2(x)')
title('Pareto front')
grid on


function y=Fun(x)
	b = [270 240 460 130];
	c = [300 300 600 200];
	t = [190 210 148 100];
	s = [200 230 160 114];
	a = [0.015 0.02 0.018 0.011];
	d = [13 13.5 14 11.5];
	y(1)=-sum(x.*(s-t));
	y(2)=sum(a.*x);
end
复制成功

最终结果展示:

最终得到多个解

x1 x2 x3 x4解得如下:

代码块
JavaScript
自动换行
复制代码
119.231391258967	0.769608488165712	0	0

119.231391258967	0.769608488165712	0	0

0.000499510291359209	120.000482867813	13.1757491344817	34.8252467588411

71.3391090591218	48.6614868846125	3.11344686170493	9.36001224382937

…… …… …… ……

27.0549008917871	92.9459248170927	9.47711506311976	25.3969178809142

8.65187243477257	111.349067997015	13.3680683558073	31.7052095195761
复制成功

        可以看到, 即使是使用智能算法, 也无法得到一般规划模型时, 我们想要得到的一个解, 最后肯定还是会出现多个解的情况. 这个时候还是需要根据多个目标函数的重要性之间的关系, 最终再进行一个取舍.

数学建模 国赛 数学模型 大学生数学建模竞赛 高教社杯 多目标规划 规划模型
本文为我原创，未经授权禁止转载
cv18490276