1. 引言
在实际生活中,各目标间的冲突性使得多目标优化问题不存在真正意义上的最优解使各子目标同时达到最优 [1] ,而只能在它们中间进行协调和折中处理,使各个子目标在给定区域都尽可能达到最优化,进而得到一个非劣解的集合。由此可见多目标问题的求解不仅仅是一个优化问题,还是一个决策问题 [2] ,当Pareto最优解集求出来之后,还需要根据决策理论,挑选出最后的折中解或最优解。
在实际应用中,解决多目标优化问题的方法主要分为两类:一种是化多为少法。化多为少的方法是将多目标问题转换为单目标问题进行求解。另外一种是分层序列法。该方法是按照一定规则将各个目标给出一个序列,每次在前一个的最优解集内求解得到下一个目标的最优解,直至求解出最终的最优解为止 [3] 。化多目标最优化问题为单目标最优化问题的一个常用方法是对各个目标函数进行线性或非线性加权,组合成一个单目标最优化问题 [4] 。
理想点法是解决多目标决策问题常用的、非常有效的求解方法。作为一种化多为少的多目标决策分析法,理想点法的基本思路是通过构造多目标问题的理想点和负理想点,并以靠近理想点和远离负理想点的程度作为决策的判断依据 [5] ,这其中需要解决的关键问题是权重向量的确定。目前,理想点法的研究方向主要分为三个方面:第一个方面是多目标决策问题中对理想点法中距离测度的研究 [6] [7] 。常见的理想点有最短距离理想点法、平方加权和理想点法和带权极大模理想点法等。第二个方面是属性权重为区间数的多目标决策问题中的理想点法的研究 [8] [9] 。这类问题主要研究各个方案与理想点的接近程度的定义方式,如理想点贴近度法、逼近理想点法、中心逼近理想点法、基于线性加权模型的理想点集成法和基于线性加权偏差平方和最小的理想点法和基于期望值的理想点法等。第三个方面是理想点法的应用方式研究,又分为理想点法的独立运用方式 [10] [11] 和融合其他多目标优化方法 [12] 的应用方式。前者主要包括多属性决策中基于理想点法的赋权方式和与多属性综合评价问题。后者主要是结合实际问题与其他的多目标优化方法进行融合 [13] [14] [15] [16] [17] 。
基于以上原因,本文提出一种包含多目标优化和多属性决策的融合相对有效性的两阶段理想点法。该方法在优化阶段采用基于相对有效性的理想点法得到一组近似的Pareto最优解,在决策阶段采用层次分析法主观赋权和变异系数法客观赋权线性组合赋权的理想点法对优化阶段获得Pareto最优解集进行排序,以辅助决策人员完成决策。
2. 多目标优化问题及求解方法
2.1. 多目标优化问题模型
多目标优化问题(VP)的标准形式为:
(1)
其中,
为约束集,也称为可行解集,r中的x称为可行解,f为目标函数,g为约束函数,x为决策变量。
多目标优化问题的最优解与单目标优化问题的最优解有着本质区别,多目标优化的本质是向量优化 [18] (Vector Programming, VP),为了正确的求解多目标优化问题,给出绝对最优解和Pareto最优解 [19] 的定义。
定义1. 设
为多目标优化模型(VP)的约束集,
是多目标优化的向量目标函数。若
,如果对所有其他的解
,
对
都成立,则称
是多目标优化模型(VP)的绝对最优解。
定义2. 设
为多目标优化模型(VP)的约束集,
是多目标优化的向量目标函数。若
,如果不存在
,使
对
都成立,且至少存在一个
使得
成立,则称
是多目标优化模型(VP)的Pareto最优解,或称Pareto有效解。
定义3. 设
为多目标优化模型(VP)的约束集,
是多目标优化的向量目标函数。若
,如果不存在
,使得对任意的
,均有
,则称
是多目标优化模型(VP)的Pareto弱有效解。
由定义可见,多目标优化问题的Pareto最优解不是唯一的。在多目标优化问题中,各目标之间冲突性使得各个子目标同时达到最优是不可能的,只能在它们中间进行协调和折中处理,使各个子目标在给定区域都尽可能达到最优化,进而得到一个非劣解的集合,这些解的特点是至少存在一个目标优于其他所有的解 [2] ,这个集合称为Pareto最优解集。
2.2. 理想点法(TOPSIS)
对于多目标问题(VP),如果决策者先能够对每个目标
给出一个值
,使其满足
(2)
则称
为理想点。特别地,如能求出
(3a)
则称
为最理想的点。理想点和最理想点有时也统称理想点。
当已知理想点
时,我们在目标空间R中,适当地引进某种模
,并考虑在这个模的意义下,目标函数
与
之间的最小“距离”,即考虑单目标问题:
(3b)
也即在(VP)的约束集合R上,寻求目标函数与理想点之间的“距离”尽可能小的解。正因为如此,我们把这种方法叫理想点法(Technique for Order Preference by Similarity to Ideal Solution, TOPSIS)。当我们给模
赋予不同的意义时,便可得到不同的理想点法。当然,模
的取法也并非随意的,一般要根据问题的实际背景或几何意义来构造。常见的理想点法有平方加权和理想点法和带权极大模理想点法等 [7] 。
2.3. 相对有效性
考虑多目标优化问题(VP),在下面的讨论中,我们将假定对任意的
,均有
。
给定(VP)的可行解
,形成如下加权形式的单目标规划决策问题:
(4)
其中
。
定义4. [20] 多目标优化决策问题(VP)的可行解
称为弱相对有效的,如果
的最优解
满足
。进一步,如果
的最优解中存在有
,并且目标值满足
,则称可行解
为相对有效的。
易知,若
为相对有效解,那么它一定也是弱相对有效解。我们称
为方案
的相对有效性评价指数。文献 [19] 中相对有效性与一般意义下的有效性概念(Pareto有效性,即Pareto最优性)之间的关系由下列定理给出:
定理1. [20] 1) 如果
为相对有效解,则
一定是Pareto有效解;2) 如果
为弱相对有效解,则
一定是弱Pareto有效解;3)如果
为弱相对有效解,且对任意的
都有
,则
一定是Pareto有效解。
定理2. [20] 如果存在Pareto有效解(弱Pareto有效解),那么一定存在相对有效解(弱相对有效解)。
定理3. [20] 如果
为下列问题的最优解:
(5)
则
为(VP)的弱相对有效解。进一步,若
,则
为(VP)的相对有效解。
3. 基于相对有效性的理想点法
理想点法是针对解的Pareto有效性设计的,它没有解决目标权重向量客观性的难题。近年来的各种新的改进的理想点法仍然没有上述问题。而相对有效性概念能够部分地克服这些缺欠。在决策过程中,计算权重和决策者参与交互进行,保证整个决策过程的客观性和决策的参与程度。
基于相对有效性 [20] 的理想点法的具体步骤如下:
第一步:确定目标函数值的取值范围及目标函数值规范化
考虑到不同目标函数下的函数值往往具有不同的量纲和量纲单位,为了消除不同量纲和量纲单位所带来的不可公度性,决策之前应对目标函数值进行无量纲化处理。于是需要依次求解,
(6)
即可作为目标
(
)的取值范围。令
(
),这时,
(
),(VP)的各目标均成为取正值的目标。易知在新的目标下,(VP)的相对有效解集不变。
第二步:由决策者给出负理想点
针对(VP),决策者对每个目标函数
给出一个值
,使其满足
(7)
则称
为负理想点。
第三步:产生供选择的解
通常情况下,各个目标函数下的函数值中的最大值不会对应同一个方案,否则,我们求解的(VP)问题已经得到解决。于是考虑到各个方案应该尽量接近理想点,尽量远离负理想点,因此,利用各目标函数值与负理想点的远离程度可建立如下的规划模型:
(8)
其中
为充分小的正向量。一般可取为
,
。
的选取是为保证所求解的相对有效性。
设最优解为
,计算
。
第四步:与决策者对话
将
和
反馈给决策者,要其作出反应,即是否对
所达到的各项目标值满意。如是,则
即为决策者满意的相对有效方案,
为对应的相对有效目标权重。否则,要决策者对下列问题作出反应:请选择你较不满意的目标值
并提出希望水平
,即希望
。把这一约束加人到原问题(VP)的规范化问题中,形成新的多目标问题并开始新的一轮迭代 [20] ,即以下列问题:
(9)
其中,
代替(VP),以下列问题:
(10)
代替(wP),转第三步。特别指出,(wP')的最优解可能会导致新的其他目标函数不能达到决策者的期望水平,所以为了避免决策者的主观意愿对决策过程的绝对控制,也为了让整个决策过程更客观合理,我们设定只有当(wP')的最优解中不满足决策者意愿的目标函数的数目小于(wP)的最优解中不满足决策者意愿的目标函数的数目时,才说明决策者给出的期望水平是合理的。否则,需要决策者对期望水平进行一定的调整。
这一过程继续,可以得到决策者满意的决策方案集,即Pareto最优解集。
4. 基于主客观线性组合赋权的理想点法
在优化阶段中,由于决策者的参与可以得到不同的Pareto最优解,进而得到一组Pareto最优解。为了使整个决策过程更合理,更多兼顾决策者的意愿和客观事实,本文对第一阶段得到的Pareto最优解集进行多属性决策,也就是对得到的Pareto最优解集进行排序。由于在多属性决策中涉及到各目标的权重,本文采用主客观线性组合赋权的方式来获得目标权重 [2] 。本文采用层次分析法进行主观赋权简单快捷,吸收了专家经验;采用变异系数法进行客观赋权具有数学基础,能够反映数据本身的信息。
设组合权向量为w,则有
。式中,
,
,
为分配给主观权向量和客观权向量的组合系数。若决策者属于经验主义者取
;若决策者属于无偏好中性者取
;若决策者属于保守主义者取
。本文采用层次分析法进行主观赋权,采用变异系数法进行客观赋权。基于主客观线性组合赋权的TOPSIS求解流程可以分为如下7步 [2] :
1) 建立规范化决策矩阵
对于第一阶段得到的m个备选方案,p个目标值的多目标决策问题,记决策方案矩阵P对应的规范化函数值矩阵为
,这里考虑的是各子目标都为效益型目标的标准形式的(VP)。进一步对初始矩
阵A规范化处理为矩阵B:
,其中
。
2) 层次分析法(Analytic Hierarchy Process, AHP)计算主观权重向量
[21] [22]
第一步:根据专家经验,对不同目标重要程度两两比较建立正互反判断矩阵 [23]
,其中
反映了第i个目标与第j个目标两两比较时相对重要性,并依据AHP理论对判断矩阵S进行一致性检验 [21] ,确定S是否具有满意一致性(即一致性比率 [24] (Consistence Rate, CR)满足
)。若满足则采用方根法计算权重向量 [25]
,反之继续到第二步。
第二步:根据下式导出依赖于矩阵S的完全一致阵 [21]
:
(11)
第三步:由完全一致阵
依据下式导出权重向量
:
(12)
3) 变异系数法计算客观权重向量
根据数据本身信息,通过下式求取各子目标的变异系数:
(13)
其中,
为是第i个目标的变异系数,也称为标准差系数;
是第i个目标的标准差;
是第i个目标的平均值。变异系数大小反映指标的离差程度和信息量,变异系数越大,赋予权重越大。变异系数权向量
。
4) 基于组合权向量w构造加权规范化矩阵C
根据
和
计算组合权向量
;利用组合权向量构造加权规范化矩阵
。
5) 确定理想点
和负理想点
理想点
,负理想点
。
6) 计算各备选方案到理想点与负理想点的欧氏距离
,
第i个备选方案到理想点的距离:
(14)
第i个备选方案到负理想点的距离:
(15)
7) 计算各方案的综合评价指数
,并按
对方案进行排序
用相对距离反映方案贴近理想点与负理想点的程度 [7]
,
越大,方案越优,反之,则方案越劣。
5. 实例演示
考虑具有四个目标函数的优化问题:
其中,
,
,
,
,
,
,
,
,
,
,
,
1) 优化阶段
步骤一:确定目标函数值的取值范围及目标函数值规范化
利用Matlab(详细代码见附录),运用线性规划和二次规划原理,得到各目标函数的取值范围分别为:
,
,
,
进而,各目标函数
规范化为
,
。
步骤二:由决策者给出负理想点
在各目标函数规范化后,假设决策者给出的负理想点为
。
步骤三:产生供选择的解
求解问题
利用Matlab,运用二次规划和非线性规划原理,得到上述问题的最优解为
,其对应的各目标函数规范化函数值为
。
步骤四:与决策者对话
假设决策者要求
,于是求解以下问题:
利用Matlab,运用二次规划和非线性规划原理,得到上述问题的最优解为
,其对应的各目标函数规范化函数值为
。与决策者的对话过程反复进行,得到决策者满意的决策方案集为矩阵P,其对应的规范化函数值矩阵为a。矩阵P,a分别为:
2) 决策阶段
步骤一:建立规范化决策矩阵
将优化阶段得到的规范化函数值矩阵a规范化处理为矩阵B:
步骤二:层次分析法计算主观权重向量
根据专家经验,对不同目标重要程度两两比较建立正互反判断矩阵
,计算得到其一致性比率
,这说明s未通过一致性检验,于是按照(11)式构造完全一致阵
,进一步由(12)式计算得到主观权重向量
。
步骤三:变异系数法计算客观权重向量
首先,以(13)式计算标准差向量和平均值向量分别为
;进而变异系数向量
;最后计算得到客观权重向量
。
步骤四:基于组合权向量w构造加权规范化矩阵C
取
并计算组合权向量
;
利用组合权向量w和B构造加权规范化矩阵
:
步骤五:确定理想点
和负理想点
计算理想点
,
负理想点
。
步骤六:计算各决策方案到理想点与负理想点的欧氏距离
,
按照(14)式和(15)式计算各决策方案到理想点
与负理想点
的欧氏距离得到距离矩阵
d的第一行为各方案与理想点
的欧氏距离,第一行为各方案与负理想点
的欧氏距离。
步骤七:计算各方案的综合评价指数
,并按
对方案进行排序
根据
计算各方案的综合评价指数得到综合评价指数向量y:
于是各方案最终排序结果为:
可以看出,方案4为最优方案,方案4对应原问题Pareto最优解
,其函数值规范化向量为
。
6. 结语
本文提出了一种融合相对有效性的两阶段理想点法。在第一阶段,运用基于相对有效性的理想点法,由于决策者的参与得到不同的Pareto最优解,决策过程中决策者参与和计算权重向量交互进行,保证了权重向量客观性和决策者的参与程度。在第二阶段,采用层次分析法主观赋权与变异系数法客观赋权线性组合赋权的理想点法对优化阶段获得的Pareto最优解排序。由于组合赋权方式吸收了主客观赋权的优点,排序结果更可靠,也更合理地反映了客观事实。
基金项目
中央高校基础研究基金“模糊信息处理与信息编码的代数理论”(2682014ZT28)。
附录
1) 优化阶段的代码
a) 求解各目标函数的取值范围的M文件BESTok如下
%% fmincon函数用来求解函数f在可行域内的最大值
A=[0,2,2,2,0,0,0,0,0,0,0;
0,-0.5,0,0,0,0,-0.25,-0.5,0,-0.25,0;
1,0,3,3,1,0,0.5,0,0,0.5,0;
0,0,0,0,0,0.25,0,0,0.25,0,0;
0,0,0,1,0,0,0,0,0.1,0.1,0;
0,0,0.5,0,0,0.5,0.5,0.5,0.4,0.3,0;
0,1,0.5,0,0,0.5,0.5,0.5,0.5,0.3,0;
0,0,0,0,1,0,0,0,0,0.3,0];
%A=[0,2,2,2,0,0,0,0,0,0,0;
%-0.5,-0.5,0,0,0,0,-0.25,-0.5,0,-0.25,0;
% 0,0,3,3,1,0,0.5,0,0,0.5,0;
%0,0,0,0,0,0.25,0,0,0.25,0,0;
% 0,0,0,1,0,0,0,0,0.1,0.1,0;
%0,0,0.5,0,0,0.5,0.5,0.5,0.4,0.3,0;
% 0,1,0.5,0,0,0.5,0.5,0.5,0.5,0.3,0;
% 1,0,0,0,1,0,0,0,0,0.3,0];
BEST=zeros(m,n);
for i=1:n
x0=A(:,i);
%x0=[0;0;1;0;0;0;0;0];
a=[-2,6,-1,1,0,0,0,0];
b=1;
aeq=[-1,-2,1,4,0,0,0,0;0,0,0,0,1,1,1,1];
beq=[1;1];
lb=[0;0;-1;1;0.1;0.1;0.1;0.1];
ub=[3;4;6;7;1;1;1;1];
%options=optimset('LargeScale','off','display','iter');
x
BEST(:,i)=x;
distancemax=-fval
end
BEST
x0
m
n
b) M文件myfun222deng0713111:
function f=myfun(x)%% x为列向量
h=[1,2,-1,1;2,5,0,3;-1,0,5,1;1,3,1,2];
f=[1;-2;0;4];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f1min=fval ;
h=-[1,2,-1,1;2,5,0,3;-1,0,5,1;1,3,1,2];
f=-[1;-2;0;4];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f1max=-fval ;
f=[-2;-1;3;-5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f2min=fval ;
f=-[-2;-1;3;-5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f2max=-fval ;
h=[1,-2,4,-1;-2,5,0,-5;4,0,7,1;-1,-5,1,2];
f=[-1;3;-7;5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f3min=fval;
h=-[1,-2,4,-1;-2,5,0,-5;4,0,7,1;-1,-5,1,2];
f=-[-1;3;-7;5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f3max=-fval ;
f=[-1;1;2;-1];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f4min=fval ;
f=-[-1;1;2;-1];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f4max=-fval;
%%函数f的定义
A1=[1,2,-1,1,0,0,0,0;2,5,0,3,0,0,0,0;-1,0,5,1,0,0,0,0;1,3,1,2,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0];
B1=[1;-2;0;4;0;0;0;0];
B2=[-2;-1;3;-5;0;0;0;0];
A3=[1,-2,4,-1,0,0,0,0;-2,5,0,-5,0,0,0,0;4,0,7,1,0,0,0,0;-1,-5,1,2,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0];
B3=[-1;3;-7;5;0;0;0;0];
B4=[-1;1;2;-1;0;0;0;0];
f=-(0.5*x'*A1*x + B1'*x-f1min)*x(5)/(f1max-f1min)-(B2'*x-f2min)*x(6)/(f2max-f2min)-(0.5*x'*A3*x + B3'*x-f3min)*x(7)/(f3max-f3min)-(B4'*x-f4min)*x(8)/(f4max-f4min);
c) M文件mycon222deng0723:
function [c,ceq]=mycon(x)%%非线性约束条件
h=[1,2,-1,1;2,5,0,3;-1,0,5,1;1,3,1,2];
f=[1;-2;0;4];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f1min=fval ;
h=-[1,2,-1,1;2,5,0,3;-1,0,5,1;1,3,1,2];
f=-[1;-2;0;4];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f1max=-fval ;
f=[-2;-1;3;-5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f2min=fval ;
f=-[-2;-1;3;-5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f2max=-fval ;
h=[1,-2,4,-1;-2,5,0,-5;4,0,7,1;-1,-5,1,2];
f=[-1;3;-7;5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f3min=fval;
h=-[1,-2,4,-1;-2,5,0,-5;4,0,7,1;-1,-5,1,2];
f=-[-1;3;-7;5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f3max=-fval ;
f=[-1;1;2;-1];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f4min=fval ;
f=-[-1;1;2;-1];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f4max=-fval;
A1=[1,2,-1,1,0,0,0,0;2,5,0,3,0,0,0,0;-1,0,5,1,0,0,0,0;1,3,1,2,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0];
B1=[1;-2;0;4;0;0;0;0];
B2=[-2;-1;3;-5;0;0;0;0];
A3=[1,-2,4,-1,0,0,0,0;-2,5,0,-5,0,0,0,0;4,0,7,1,0,0,0,0;-1,-5,1,2,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0];
B3=[-1;3;-7;5;0;0;0;0];
B4=[-1;1;2;-1;0;0;0;0];
c(1)=f1min-(0.5*x'*A1*x + B1'*x);
c(2)=0.5*x'*A1*x + B1'*x-f1max;
c(3)=f2min-B2'*x;
c(4)=B2'*x-f2max;
c(5)=f3min-(0.5*x'*A3*x + B3'*x);
c(6)=0.5*x'*A3*x + B3'*x-f3max;
c(7)=f4min-B4'*x;
c(8)=B4'*x-f4max;
%c(9)=-(0.5*x'*A1*x + B1'*x-f1min)/(f1max-f1min)+0.8;
%c(10)=-(B2'*x-f2min)/(f2max-f2min)+0.6;
%c(11)=-(0.5*x'*A3*x + B3'*x-f3min)/(f3max-f3min)-0.8;
%c(12)=-(B4'*x-f4min)/(f4max-f4min)+0.6;
ceq=[];
%c=[];
d) 规范化最优解矩阵M文件Percentok
h=[1,2,-1,1;2,5,0,3;-1,0,5,1;1,3,1,2];
f=[1;-2;0;4];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f1min=fval
h=-[1,2,-1,1;2,5,0,3;-1,0,5,1;1,3,1,2];
f=-[1;-2;0;4];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f1max=-fval
f=[-2;-1;3;-5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f2min=fval
f=-[-2;-1;3;-5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f2max=-fval
h=[1,-2,4,-1;-2,5,0,-5;4,0,7,1;-1,-5,1,2];
f=[-1;3;-7;5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f3min=fval
h=-[1,-2,4,-1;-2,5,0,-5;4,0,7,1;-1,-5,1,2];
f=-[-1;3;-7;5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f3max=-fval
f=[-1;1;2;-1];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f4min=fval
f=-[-1;1;2;-1];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f4max=-fval
A1=[1,2,-1,1,0,0,0,0;2,5,0,3,0,0,0,0;-1,0,5,1,0,0,0,0;1,3,1,2,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0];
B1=[1;-2;0;4;0;0;0;0];
B2=[-2;-1;3;-5;0;0;0;0];
A3=[1,-2,4,-1,0,0,0,0;-2,5,0,-5,0,0,0,0;4,0,7,1,0,0,0,0;-1,-5,1,2,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0];
B3=[-1;3;-7;5;0;0;0;0];
B4=[-1;1;2;-1;0;0;0;0];
BEST=[2.21092.59562.88932.78302.36512.64962.21082.64502.79942.21173.0000;
0.90810.99621.06801.04160.94251.00920.90811.00811.04570.90831.0958;
1.02710.94670.90850.92030.99010.93831.02710.93900.91841.02690.8982;
1.00001.16031.27921.23651.06501.18241.00001.18051.24311.00041.3234;
0.10000.10000.10000.10000.10000.10000.10000.10000.10000.10000.1000;
0.10000.10000.10000.10000.10000.10000.10000.10000.10000.10000.1000;
0.70000.70000.70000.70000.70000.70000.70000.70000.70000.70000.7000;
0.10000.10000.10000.10000.10000.10000.10000.10000.10000.10000.1000]
percent=zeros(n,4);%percent为最优解对应原各目标函数的最优值的归一化矩阵
for i=1:n
x=BEST(:,i);
percent(i,1)=(0.5*x'*A1*x + B1'*x-f1min)/(f1max-f1min);
percent(i,2)=(B2'*x-f2min)/(f2max-f2min);
percent(i,3)=(0.5*x'*A3*x + B3'*x-f3min)/(f3max-f3min);
percent(i,4)=(B4'*x-f4min)/(f4max-f4min);
end
percent
2) 决策阶段的代码
a) 求解最后决策方案的M文件trydengokjieguook0815
disp('=====================================================');
%BEST =
% Columns 1 through 9
% 2.2109 2.5956 2.8893 2.7830 2.3651 2.6496 2.2108 2.6450 2.7994
% 0.9081 0.9962 1.0680 1.0416 0.9425 1.0092 0.9081 1.0081 1.0457
% 1.0271 0.9467 0.9085 0.9203 0.9901 0.9383 1.0271 0.9390 0.9184
% 1.0000 1.1603 1.2792 1.2365 1.0650 1.1824 1.0000 1.1805 1.2431
% 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000
% 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000
% 0.7000 0.7000 0.7000 0.7000 0.7000 0.7000 0.7000 0.7000 0.7000
% 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000
% Columns 10 through 11
% 2.2117 3.0000
% 0.9083 1.0958
% 1.0269 0.8982
% 1.0004 1.3234
% 0.1000 0.1000
% 0.1000 0.1000
% 0.7000 0.7000
% 0.1000 0.1000
disp('第一阶段优化得到的最优解矩阵:');%BEST为第一阶段(融合相对有效性的理想点法)的一系列最优解(运行 BESTok.m 得到)
BEST = [2.21092.59562.88932.78302.36512.64962.21082.64502.79942.21173.0000;
0.90810.99621.06801.04160.94251.00920.90811.00811.04570.90831.0958;
1.02710.94670.90850.92030.99010.93831.02710.93900.91841.02690.8982;
1.00001.16031.27921.23651.06501.18241.00001.18051.24311.00041.3234;
0.10000.10000.10000.10000.10000.10000.10000.10000.10000.10000.1000;
0.10000.10000.10000.10000.10000.10000.10000.10000.10000.10000.1000;
0.70000.70000.70000.70000.70000.70000.70000.70000.70000.70000.7000;
0.10000.10000.10000.10000.10000.10000.10000.10000.10000.10000.1000]
disp('最优解对应原各目标函数的最优值的归一化矩阵:');%percent为最优解对应原各目标函数的最优值的归一化矩阵 (运行 将BEST代入Percentok.m 得到)
percent =[0.3687 0.7380 1.0000 0.6155;
0.6009 0.6052 1.0000 0.5522;
0.8035 0.5096 1.0000 0.5094;
0.7277 0.5437 1.0000 0.5244;
0.4571 0.6836 1.0000 0.5890;
0.6366 0.5872 0.9999 0.5440;
0.3687 0.7380 1.0000 0.6155;
0.6335 0.5888 1.0000 0.5447;
0.7393 0.5384 1.0000 0.5221;
0.3692 0.7377 1.0000 0.6154;
0.8851 0.4745 1.0000 0.4942]
%最小最大值规范化矩阵percent
maxpercent=max(percent)
minpercent=-max(-percent)
MM=zeros(M,N);
for i=1:M
for j=1:N
MM(i,j)=(percent(i,j)-minpercent(1,j))/(maxpercent(1,j)-minpercent(1,j));
end
end
disp('规范化函数值矩阵percent的规范化矩阵为:');
MM
disp('层次分析法得到的主观权重向量为:');%(运行 Ahpweight2220815ok.m 得到)
weight1 =[0.4821 0.2913 0.1158 0.1107]
disp('变异系数法得到的客观权重向量为:');%(运行 将MM代入bianyixishufa222ok.m 得到)
weight2 =[0.3151 0.2738 0.1303 0.2808]
weight=zeros(1,4);
for i=1:4
weight(1,i)=0.5*weight1(1,i)+0.5*weight2(1,i);
end
disp('加权组合权重向量为:');
weight
disp('加权规范化矩阵为:');
weightpercent=zeros(m,n);
for i=1:m
for j=1:n
weightpercent(i,j)=percent(i,j)*weight(1,j);
end
end
weightpercent
disp('正理想点为:');
posideal=max(weightpercent)
disp('负理想点为:');
p=-weightpercent;
negideal=-max(p)
disp('各方案到正负理想点的距离矩阵为:');
distance=zeros(m,3); %distance为各方案到正负理想点的距离及相对距离矩阵
sum1=0;
sum2=0;
for i=1:m
for j=1:n
sum1=sum1+(weightpercent(i,j)-posideal(1,j))*(weightpercent(i,j)-posideal(1,j));
sum2=sum2+(weightpercent(i,j)-negideal(1,j))*(weightpercent(i,j)-negideal(1,j));
end
distance(i,1)=sqrt(sum1);
distance(i,2)=sqrt(sum2);
distance(i,3)=distance(i,2)/distance(i,1);
end
distance
disp('相对距离排序');
aa=Y'
disp('对应的方案排序');
bb=I'
fprintf('\n方案%i为最优决策方案\n',bb(1,n));
h=[1,2,-1,1;2,5,0,3;-1,0,5,1;1,3,1,2];
f=[1;-2;0;4];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f1min=fval ;
h=-[1,2,-1,1;2,5,0,3;-1,0,5,1;1,3,1,2];
f=-[1;-2;0;4];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f1max=-fval ;
f=[-2;-1;3;-5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f2min=fval ;
f=-[-2;-1;3;-5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f2max=-fval ;
h=[1,-2,4,-1;-2,5,0,-5;4,0,7,1;-1,-5,1,2];
f=[-1;3;-7;5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f3min=fval;
h=-[1,-2,4,-1;-2,5,0,-5;4,0,7,1;-1,-5,1,2];
f=-[-1;3;-7;5];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f3max=-fval ;
f=[-1;1;2;-1];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f4min=fval ;
f=-[-1;1;2;-1];
a=[-2,6,-1,1];
b=1;
aeq=[-1,-2,1,4];
beq=1;
lb=[0;0;-1;1];
ub=[3;4;6;7];
f4max=-fval;
A1=[1,2,-1,1,0,0,0,0;2,5,0,3,0,0,0,0;-1,0,5,1,0,0,0,0;1,3,1,2,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0];
B1=[1;-2;0;4;0;0;0;0];
B2=[-2;-1;3;-5;0;0;0;0];
A3=[1,-2,4,-1,0,0,0,0;-2,5,0,-5,0,0,0,0;4,0,7,1,0,0,0,0;-1,-5,1,2,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0];
B3=[-1;3;-7;5;0;0;0;0];
B4=[-1;1;2;-1;0;0;0;0];
BEST=[2.21092.59562.88932.78302.36512.64962.21082.64502.79942.21173.0000;
0.90810.99621.06801.04160.94251.00920.90811.00811.04570.90831.0958;
1.02710.94670.90850.92030.99010.93831.02710.93900.91841.02690.8982;
1.00001.16031.27921.23651.06501.18241.00001.18051.24311.00041.3234;
0.10000.10000.10000.10000.10000.10000.10000.10000.10000.10000.1000;
0.10000.10000.10000.10000.10000.10000.10000.10000.10000.10000.1000;
0.70000.70000.70000.70000.70000.70000.70000.70000.70000.70000.7000;
0.10000.10000.10000.10000.10000.10000.10000.10000.10000.10000.1000];
fprintf('\n方案%i的最优值及其归一化值分别为:\n',bb(1,n));
bestok=BEST(:,bb(1,n))
percentok=zeros(1,4);%percent为最优解对应原各目标函数的最优值的归一化矩阵
x=BEST(:,bb(1,n));
percentok(1,1)=(0.5*x'*A1*x + B1'*x-f1min)/(f1max-f1min);
percentok(1,2)=(B2'*x-f2min)/(f2max-f2min);
percentok(1,3)=(0.5*x'*A3*x + B3'*x-f3min)/(f3max-f3min);
percentok(1,4)=(B4'*x-f4min)/(f4max-f4min);
percentok
b) 计算主观权重向量M文件Ahpweight2220815ok
disp('========================================');
%A=[1 2 6;1/2 1 2;1/6 1/2 1];
A=[1,3,5,2;1/3,1,6,2;1/5,1/6,1,3;1/2,1/2,1/3,1];
disp('DecMatrix=A');
%disp('一致性比率为:');
%CR
if CR<0.1
disp('判断矩阵A通过一致性检验');
M=ones(1,m0);
M1=ones(1,m0);
weight=ones(1,m0);
for i=1:m0 %%计算判断矩阵每一行元素乘积M(1,i)
for j=1:m0
M(1,i)=M(1,i)*A(i,j);
end
end
for i=1:m0
%M1(1,i)=nthroot(M(1,i), m0);
M1(1,i)=power(M(1,i), 1/m0);
end
sum=0;
for i=1:m0
sum=sum+M1(1,i);
end
for i=1:m0
weight(1,i)= M1(1,i)/sum;
end
disp('采用方根法得到的权重向量为:');
weight
else
disp('判断矩阵A未通过一致性检验');
A1=zeros(m0,m0);%%修正判断矩阵
for i=1:m0
for j=1:m0
m=1;
for k=1:m0
m=m*A(i,k)*A(k,j);
end
A1(i,j)=power(m, 1/m0);
%A1(i,j)=nthroot(m, m0);
end
end
disp('将判断矩阵A修正为A1:');
A1
sum1=zeros(1,m0);
weight1=zeros(1,m0);
for i=1:m0
for j=1:m0
sum1(1,i)=sum1(1,i)+A1(j,i);
end
weight1(1,i)=1/sum1(1,i);
end
disp('权重向量为:');
weight1
end
c) 计算客观权重向量M文件bianyixishufa222ok
disp('================================================');
%A=[1 2 3;1 1 1]
%这里A为第一阶段得到的规范化函数值矩阵percent的规范化矩阵(最优解对应原各目标函数的最优值归一化矩阵)
percent =[0.3687 0.7380 1.0000 0.6155;
0.6009 0.6052 1.0000 0.5522;
0.8035 0.5096 1.0000 0.5094;
0.7277 0.5437 1.0000 0.5244;
0.4571 0.6836 1.0000 0.5890;
0.6366 0.5872 0.9999 0.5440;
0.3687 0.7380 1.0000 0.6155;
0.6335 0.5888 1.0000 0.5447;
0.7393 0.5384 1.0000 0.5221;
0.3692 0.7377 1.0000 0.6154;
0.8851 0.4745 1.0000 0.4942];
A=[ 0 1.0000 1.0000 1.0000;
0.4497 0.4960 1.0000 0.4782;
0.8420 0.1332 1.0000 0.1253;
0.6952 0.2626 1.0000 0.2490;
0.1712 0.7935 1.0000 0.7815;
0.5188 0.4277 0 0.4106;
0 1.0000 1.0000 1.0000;
0.5128 0.4338 1.0000 0.4163;
0.7177 0.2425 1.0000 0.2300;
0.0010 0.9989 1.0000 0.9992;
1.0000 0 1.0000 0];
%standard=std(A,0)%除以n-1 样本标准差
standard=std(A,1)%除以n 标总体准差
aver=zeros(1,n1);
stdnumber=zeros(1,n1);
for i=1:n1
aver(1,i)=mean(A(:,i));
stdnumber(1,i)=standard(1,i)/mean(A(:,i));
%stdnumber(1,i)=standard(1,i)/aver(1,i);
end
aver
stdnumber
sum=0;
for i=1:n1
sum=sum+stdnumber(1,i);
end
weight2=zeros(1,n1);
for i=1:n1
weight2(1,i)=stdnumber(1,i)/sum;
end
disp('客观权重向量为:');
weight2