參考資料

魯棒最佳化（Robust optimization）

魯棒最佳化是指，定義問題的資料不正確或者不確定的情況下，可以返回信賴的結果的最佳化問題的

建模

技術及演算法。

由於最佳化的結果對於外部的擾亂表現的非常不穩定，存在完全無法信賴的結果的情況。比如（Aharon Ben-Tal and Arkadi Nemirovski ），非常有名的

線性規劃

問題的 benchmark （NETLIB）中，對定義的線性約束的係數做比較小的擾動的情況下， 問題可能變得不可行。

二次錐最佳化問題（Second-order cone optimization problem）

這類最佳化問題，如他的名字一樣，因為約束裡有二次錐約束。

最基本的3次元的二次錐為

\{(z,x,y):z\geq \sqrt{x^{2}+y^{2}}\}

魯棒線性規劃問題的建模及求解

3次元的二次錐

看它長得像什麼呢？冰激凌，所以叫它冰激凌錐。別急它還有個更加高大上的名字，叫做洛倫茲錐(Lorentz cone)。錐的種類很多，大家可詳細搜尋一下。

混合問題（Product mix problem）的建模

這裡我們考慮一個古典的線性規劃問題-混合問題，它可描述為：

這裡有一個工廠。它要把四種原料混合之後，製造一種產品。對於原料包含三種成分。我們想混合成為，使成分1百分之20以上，成分2百分之30以上，成分3百分之20以上。各原料成分的含有比例如圖，原料單價每噸分別是5，6，8，20元。那麼，怎麼混合原料使得用最小費用來製造產品？

魯棒線性規劃問題的建模及求解

各原料成分含有比例

我們可把混合問題用線性規劃問題進行建模如下：

\begin{eqnarray} (\clubsuit) \left| \begin{array}{ll} \mbox{min} & \sum_{i=1}^{4} p_{i}x_{i}\\ \mbox{s.t} & \sum_{i=1}^{4} x_{i}=1\\ &\sum_{i=1}^{4}a_{i,k}x_{i}\geq LB_{k}\quad k=1,2,3\\ & x_{i} \geq 0\quad i=1,2,3,4 \nonumber \end{array} \right. \end{eqnarray}

(\clubsuit)

中，原料

i

的價格是

p_{i}

，在原料

i

裡包含成分

k

的比率為

a_{i,k}

，產品中應該包含的成分

k

的比率的下限為

LB_{k}

。表示原料

i

的混合比率的變數為

x_{i}

。這樣用求解器直接求解即可。

但是，

實際中，原料裡包含的成分的比率能夠正確的確定的情況很少，多少會伴隨有點誤差。

在這裡我們考慮原料

i

裡包含成分

k

的比率

a_{i,k}

伴隨有誤差

e_{i,k}

。另外，我們需要對這個誤差在一定範圍內定義。這裡我們考慮

\sum_{i=1}^{4}e_{i,k}^{2}\leq \epsilon^{2}

\epsilon

表示限制誤差量的引數。

那麼，考慮上述的誤差，

(\clubsuit)

中的約束條件

\sum_{i=1}^{4}a_{i,k}x_{i}\geq LB_{k}

可重新整理為

\sum_{i=1}^{4}(a_{i,k}+e_{i,k})x_{i}\geq LB_{k}

進而可得

\min \{\sum_{i=1}^{4}(a_{i,k}+e_{i,k})x_{i}:\sum_{i=1}^{4}e_{i,k}^{2}\leq \epsilon^{2}\}\geq LB_{k}

或

\min \{\sum_{i=1}^{4}e_{i,k}x_{i}:\sum_{i=1}^{4}e_{i,k}^{2}\leq \epsilon^{2}\}\geq LB_{k}-\sum_{i=1}^{4}a_{i,k}x_{i}

這裡，在式子左邊，取

\min

的地方是表示在4次元的半徑為

\epsilon

的球上的

線性函式

\sum_{i=1}^{4}e_{i,k}x_{i}

的最小化。但是，變數是

e_{i,k}

。

什麼？什麼？4次元!

為了問題的簡單化，我們考慮2次元的半徑為1的球（圓）上的線性函式

e_{1}v_{1}+e_{2}v_{2}

最小化的情況。

魯棒線性規劃問題的建模及求解

從上圖可知，在圓上使

e_{1}v_{1}+e_{2}v_{2}

最小的點

(e_{1},e_{2})

有

(e_{1},e_{2})=-1\times\frac{(v_{1},v_{2})}{\sqrt{v_{1}^{2}+v_{2}^{2}}}

進而推廣到現在這個問題上，可得線性函式在球上的最優解為

e_{i,k}=-\epsilon\times\frac{x_{i}}{\|x \|} \quad i=1,2,3,4

但是，

\|x\|=\sqrt{\sum_{i=1}^{4}x_{i}^{2}}

。所以，我們的約束進而可表示為

\min \{\sum_{i=1}^{4}e_{i,k}x_{i}:\sum_{i=1}^{4}e_{i,k}^{2}\leq \epsilon^{2}\}=-\epsilon\sum_{i=1}^{4}\frac{x_{i}^{2}}{\|x\|}=-\epsilon \|x\|

最終，

(\clubsuit)

中含有誤差的約束可表示為二次錐約束，如下

\epsilon \|x\| \leq -LB_{k}+\sum_{i=1}^{4}a_{i,k}x_{i}

Reformulation後，魯棒最佳化問題可用二次錐最最佳化問題建模為

\begin{eqnarray} (\spadesuit) \left| \begin{array}{ll} \mbox{min} & \sum_{i=1}^{4} p_{i}x_{i}\\ \mbox{s.t} & \sum_{i=1}^{4} x_{i}=1\\ &\sqrt{\epsilon^{2}\sum_{i=1}^{4}x_{i}^{2}}\leq -LB_{k}+\sum_{i=1}^{4}a_{i,k}x_{i}\quad k=1,2,3\\ & x_{i} \geq 0\quad i=1,2,3,4 \nonumber \end{array} \right. \end{eqnarray}

如果這裡

\epsilon=0

那麼

(\spadesuit)=(\clubsuit)

。

(\spadesuit)

的求解

(\spadesuit)

為二次錐最佳化問題，Gurobi （內點法已經OK了，From Version 5。0）是可以 handle。這裡，我們用程式語言Python呼叫Gurobi的API進行求解。

#code reference： Prof。 Kubo

from

gurobipy

import

*

def

prodmix

（

I

，

K

，

a

，

p

，

epsilon

，

LB

）：

“”“prodmix： robust production planning using soco

Parameters：

I - set of materials

K - set of components

a［i］［k］ - coef。 matrix

p［i］ - price of material i

LB［k］ - amount needed for k

Returns a model， ready to be solved。

”“”

model

=

Model

（

“robust product mix”

）

x

，

rhs

=

{}，{}

for

i

in

I

：

x

［

i

］

=

model

。

addVar

（

vtype

=

“C”

，

name

=

“x（

%s

）”

%

i

）

for

k

in

K

：

rhs

［

k

］

=

model

。

addVar

（

vtype

=

“C”

，

name

=

“rhs（

%s

）”

%

k

）

model

。

update

（）

model

。

addConstr

（

quicksum

（

x

［

i

］

for

i

in

I

）

==

1

）

for

k

in

K

：

model

。

addConstr

（

rhs

［

k

］

==

-

LB

［

k

］

+

quicksum

（

a

［

i

，

k

］

*

x

［

i

］

for

i

in

I

）

）

model

。

addConstr

（

quicksum

（

epsilon

*

epsilon

*

x

［

i

］

*

x

［

i

］

for

i

in

I

）

<=

rhs

［

k

］

*

rhs

［

k

］）

model

。

setObjective

（

quicksum

（

p

［

i

］

*

x

［

i

］

for

i

in

I

），

GRB

。

MINIMIZE

）

model

。

update

（）

model

。

__data

=

x

，

rhs

return

model

def

make_data

（）：

a

=

{

（

1

，

1

）：

。

25

，

（

1

，

2

）：

。

15

，

（

1

，

3

）：

。

2

，

（

2

，

1

）：

。

3

，

（

2

，

2

）：

。

3

，

（

2

，

3

）：

。

1

，

（

3

，

1

）：

。

15

，

（

3

，

2

）：

。

65

，

（

3

，

3

）：

。

05

，

（

4

，

1

）：

。

1

，

（

4

，

2

）：

。

05

，

（

4

，

3

）：

。

8

}

epsilon

=

0。01

#0。02

I

，

p

=

multidict

（{

1

：

5

，

2

：

6

，

3

：

8

，

4

：

20

}）

K

，

LB

=

multidict

（{

1

：

。

2

，

2

：

。

3

，

3

：

。

2

}）

return

I

，

K

，

a

，

p

，

epsilon

，

LB

I

，

K

，

a

，

p

，

epsilon

，

LB

=

make_data

（）

model

=

prodmix

（

I

，

K

，

a

，

p

，

epsilon

，

LB

）

model

。

optimize

（）

print

（

“obj：”

，

model

。

ObjVal

）

x

，

rhs

=

model

。

__data

for

i

in

I

：

print

（

i

，

x

［

i

］

。

X

）

如果感興趣的可以做一下當

\epsilon>0

，

(\spadesuit)

和

(\clubsuit)

的最優解以及

目標函式

值的差別。

如果有什麼疑問或者建議，可以給我發郵件。➡ zhaoyou728 atoutlook。com