author:旭寶ww
DateTime:2020/7/2
一、引言
對于多于一個的目标函數在給定區域上的最優化問題稱為多目标規劃問題。在多目标規劃中,各目标之間是互相沖突的,不一定存在所有目标上都是最優的解。是以多目标問題的解構成一個集合,他們之間不能簡單地比較好壞,這樣的解稱為非支配解(有效解) 或者 Pareto最優解。
注意:多目标規劃不同于單目标規劃,在數學模組化的結果中不應當給出一個最優解,Pareto最優解應當是一組解,以供甲方選擇。
二、模型執行個體
1 一般地,多目标規劃的數學模型如下:
m i n F ( x ) = ( f 1 ( x ) , f 2 ( x ) , … , f m ( x ) ) s . t . x ∈ Ω \mathbf{min \quad F(x)=\left (f_1(x),f_2(x),…,f_m(x)\right) }\\ \mathbf{s.t.\quad x\in \Omega} minF(x)=(f1(x),f2(x),…,fm(x))s.t.x∈Ω
其中 x = ( x 1 , x 2 , … , x n ) x=(x_1,x_2,…,x_n) x=(x1,x2,…,xn)所在的空間 Ω \Omega Ω稱為問題的決策空間(可行解空間),m維向量 F ( x ) F(x) F(x)所在的空間稱為問題的目标空間。
其中 Ω = { x ∈ R n ∣ g i ( x ) ≤ 0 , i = 1 , 2 , … , p } \Omega =\left \{x\in R^{n} |g_i(x)\le 0,i=1,2,…,p \right \} Ω={x∈Rn∣gi(x)≤0,i=1,2,…,p}
2 相關概念
定義1:對最小化問題,一個向量 u = ( u 1 , u 2 , … , u m ) u=(u_1,u_2,…,u_m) u=(u1,u2,…,um)支配另一個向量 v = ( v 1 , v 2 , … , v m ) v=(v_1,v_2,…,v_m) v=(v1,v2,…,vm),當且僅當
{ u i ≤ v i , ∀ i ∈ { 1 , 2 , … , m } ∃ j ∈ { 1 , 2 , … , m } \begin{cases} u_i \le v_i,\forall i \in \left \{ 1,2,…,m \right \} \\ \exists j \in \left \{ 1,2,…,m \right \} \end{cases} {ui≤vi,∀i∈{1,2,…,m}∃j∈{1,2,…,m}
定義2:對于任意兩個自變量向量 x 1 , x 2 ∈ Ω x_1,x_2\in \Omega x1,x2∈Ω,如果向量 ( f 1 ( x 1 ) , f 2 ( x 1 ) , … , f m ( x 1 ) ) (f_1(x_1),f_2(x_1),…,f_m(x_1)) (f1(x1),f2(x1),…,fm(x1))支配向量 ( f 1 ( x 2 ) , f 2 ( x 2 ) , … , f m ( x 2 ) ) (f_1(x_2),f_2(x_2),…,f_m(x_2)) (f1(x2),f2(x2),…,fm(x2)),則稱 x 1 x_1 x1支配 x 2 x_2 x2。
定義3:如果 Ω \Omega Ω中沒有支配 x x x的解,則稱 x x x是問題的一個Pareto最優解,Pareto最優解的集合稱為Pareto最優解集,Pareto最優解集在目标空間的像集稱為Pareto前沿(Pareto Front)。
例如,此圖的橫坐标為 F 1 ( x ) F_1(x) F1(x),縱坐标為 F 2 ( x ) F_2(x) F2(x),紅色的連線即為該問題的Pareto前沿,H和G兩點的自變量 x x x都是此問題的解。
三、求解方法
多目标優化問題的算法大體分為傳統的數學規劃方法和智能優化方法。傳統方法隻能得到非支配解,故需要重點掌握智能優化方法的MATLAB實作。
3.1傳統的數學規劃方法
1.主要目标法:确定一個目标為主要目标,确定合适的界限值把次要目标作為限制條件。
{ m i n f 1 ( x ) s . t . f i ( x ) ≤ a i , i = 2 , 3 , … , m x ∈ Ω \begin{cases} min \quad f_1(x)\\ s.t. \quad f_i(x) \le a_i,i=2,3,…,m\\ \quad\quad x \in \Omega \end{cases} ⎩⎪⎨⎪⎧minf1(x)s.t.fi(x)≤ai,i=2,3,…,mx∈Ω
其中界限值 a i ≥ min x ∈ Ω f i ( x ) , i = 2 , 3 , … , m a_i\ge \min \limits_{x\in \Omega} f_i(x),i=2,3,…,m ai≥x∈Ωminfi(x),i=2,3,…,m
2.分層序列法:先對m個目标的重要性排序。先求問題:
P ( 1 ) { m i n f 1 ( x ) s . t . x ∈ Ω P(1)\quad \begin{cases} min \quad f_1(x)\\ s.t.\quad x \in \Omega \end{cases} P(1){minf1(x)s.t.x∈Ω
的最優解 x ( 1 ) x^{(1)} x(1)和最優值 f 1 ∗ f_1^{*} f1∗。
再求問題:
P ( 2 ) { m i n f 2 ( x ) s . t . x ∈ Ω 1 = Ω ⋂ { x ∣ f 1 ( x ) ≤ f 1 ∗ } P(2)\quad \begin{cases} min \quad f_2(x)\\ s.t.\quad x \in \Omega_1=\Omega \bigcap \left \{ x|f_1(x) \le f_1^* \right \} \end{cases} P(2){minf2(x)s.t.x∈Ω1=Ω⋂{x∣f1(x)≤f1∗}
的最優解 x ( 2 ) x^{(2)} x(2)和最優值 f 2 ∗ f_2^{*} f2∗。
如此進行下去,直到求出第 m m m個問題:
P ( m ) { m i n f m ( x ) s . t . x ∈ Ω m − 1 = Ω m − 2 ⋂ { x ∣ f m − 1 ( x ) ≤ f m − 1 ∗ } P(m)\quad \begin{cases} min \quad f_m(x)\\ s.t.\quad x \in \Omega_{m-1}=\Omega_{m-2} \bigcap \left \{ x|f_{m-1}(x) \le f_{m-1}^* \right \} \end{cases} P(m){minfm(x)s.t.x∈Ωm−1=Ωm−2⋂{x∣fm−1(x)≤fm−1∗}
的最優解 x ( m ) x^{(m)} x(m)和最優值 f m ∗ f_m^{*} fm∗。
于是 x ∗ = x ( m ) x^*=x^{(m)} x∗=x(m)就是多目标規劃問題的一個非支配解。
3.權重法:對m個目标加以适當的權重,化作單目标規劃
{ m i n ∑ i = 1 m ω i f i ( x ) s . t . x ∈ Ω \begin{cases} min \sum\limits_{i=1}^{m}\omega_if_i(x)\\ s.t.\quad x \in \Omega \end{cases} ⎩⎨⎧mini=1∑mωifi(x)s.t.x∈Ω
4.MATLAB求解單目标規劃問題
[x,fval]=quadprog(H,f,A,b,Aeq,beq,lb,ub) %%求解二次規劃 x T A x x^TAx xTAx
[x,fval]=fmincon(f,x0,A,b,Aeq,beq,lb,ub,nonlcon) %%求解限制下的多元函數最小值
3.2 智能優化算法(NSGA-Ⅱ)
能力有限,無法了解算法的原理,貼一位大佬的部落格:
https://blog.csdn.net/qq_40434430/article/details/82876572
MATLAB中的内置函數gamultiobj() (ga:遺傳;multiobj:多目标 )是基于NSGA-Ⅱ改進的一種多目标遺傳算法。
gamultiobj适用于求解以下形式的多目标問題:
{ m i n F ( x ) s . t . A X ≤ b A e q ∗ X = b e q l b ≤ X ≤ u b \begin{cases} min \quad F(x)\\ s.t.\quad AX\le b\\ \qquad Aeq*X=beq\\ \qquad lb\le X\le ub \end{cases} ⎩⎪⎪⎪⎨⎪⎪⎪⎧minF(x)s.t.AX≤bAeq∗X=beqlb≤X≤ub
調用格式:
[X,FVAL]=gamultiobj( f, NVARS, A, b, Aeq, beq, lb, ub, nonlcon, options)
%X是Pareto解集
%FVAL是目标函數在Pareto解上的函數值
%f是需要計算的目标函數
%NVARS是變量的個數
%options是繪制Pareto前沿的函數gaoptimset的參數
執行個體:
{ m i n f 1 ( x ) = x 1 4 − 10 x 1 2 + x 1 x 2 + x 2 4 − x 1 2 x 2 2 m i n f 2 ( x ) = x 1 4 + x 1 x 2 + x 2 4 − x 1 2 x 2 2 s . t . − 5 ≤ x 1 , x 2 ≤ 5 \begin{cases} min\quad f_1(x)=x_1^4-10x_1^2+x_1x_2+x_2^4-x_1^2x_2^2\\ min\quad f_2(x)=x_1^4+x_1x_2+x_2^4-x_1^2x_2^2\\ s.t.\quad -5\le x_1,x_2\le 5 \end{cases} ⎩⎪⎨⎪⎧minf1(x)=x14−10x12+x1x2+x24−x12x22minf2(x)=x14+x1x2+x24−x12x22s.t.−5≤x1,x2≤5
function y=Fun(x)
y(1)=x(1)^4-10*x(1)^2+x(1)*x(2)+x(2)^4-x(1)^2*x(2)^2;
y(2)=x(2)^4-x(1)^2*x(2)^2+x(1)^4+x(1)*x(2);
end
clear;clc;
fitnessfcn=@Fun %%調用函數句柄
nvars=2; %%變量個數
lb=[-5,-5];
ub=[5,5];
A=[];b=[];
Aeq=[];beq=[];
options=gaoptimset('paretoFraction',0.3,'populationsize',100,'generations',200,'stallGenLimit',200,'TolFun',1e-10,'PlotFcns',@gaplotpareto);
%最優個體系數paretoFraction為0.3
%種群大小populationsize為100
%最大進化代數generations為200
%停止代數stallGenLimit為200
%适應度函數偏差TolFun為1e-10
%函數gaplotpareto:繪制Pareto前沿
[x,fval]=gamultiobj(fitnessfcn,nvars,A,b,Aeq,beq,lb,ub,options)
MATLAB運作結果如下: