floptで双対問題の作成とその求解

本稿ではfloptを用いて双対問題の作成のその求解を行います。

双対問題
実装
線形計画法の双対問題

双対問題

min  f(x)
s.t. g_i(x) <= 0  (for all i)
     x is continuous vector

このような問題の場合、ラグランジュ関数を用いた双対問題を次のように構築することができます。まずラグランジュ関数を

$L(\theta, x) = f(x) + \sum_i \theta_i g_i(x)$

と定義します。ここで $\theta$ は非負の値を取る変数です。このラグランジュ関数を用いて

$max_\theta min_x L(\theta, x)$

という問題が定義できるのですが、これが主問題の双対問題になっています。すなわち主問題の最適値を $P^*$ とすると $P^* \geq \min_x L(\theta, x) (\forall \theta \in \mathbb{R}_{+})$ が成立します。従って、この双対問題を用いることで、主問題の最適値の下限を求めることができます。

実装

min  (x-1)^2 + (y-2.5)^2
s.t. x - 2y >= -2

この主問題を用いて実装の確認を行います。この問題の最適解は(x, y) = (1.4, 1.7)で最適値は0.8です。

import flopt

x = flopt.Variable("x", lowBound=0)
y = flopt.Variable("y", lowBound=0)

primary_prob = flopt.Problem(name="Primary problem", sense="Minimize")
primary_prob += (x - 1) ** 2 + (y - 2.5) ** 2
primary_prob += x - 2 * y >= -2
primary_prob.show()
>>> Name: Primary problem
>>>   Type         : Problem
>>>   sense        : Minimize
>>>   objective    : (x-1)^2+(y-2.5)^2
>>>   #constraints : 1
>>>   #variables   : 2 (Continuous 2)
>>> 
>>>   C 0, name None, -2-(x-(2*y)) <= 0
>>> 
>>>   V 0, name x, Continuous 0 <= x <= None
>>>   V 1, name y, Continuous 0 <= y <= None

実際にfloptで解いてみると、

primary_prob.solve()
print("primary obj", primary_prob.getObjectiveValue())
for var in primary_prob.getVariables():
    print(var.name, var.value())
>>> primary obj 0.8000021696245797
>>> x 1.4000008269334314
>>> y 1.699999057452336

と、最適解が得られていますね。

制約を不等式に変換

Problem::boundsToIneqを用いて変数の制約をすべて不等式制約に変換します。すなわち、a <= x <= b という変数xの範囲を、a - x <= 0 と x - b <= 0の二つの制約に変換します。また、同時にここで最小化問題に変換しておきましょう。

prob = primary_prob.clone(variable_clone=True)  # copy problem
prob = prob.boundsToIneq()  # convert bounds of variables to constraints
if prob.sense.lower() == "maximize":
    prob.obj = -prob.obj
    prob.sense = "Minimize"
prob.show()
>>> Name: Primary problem
>>>   Type         : Problem
>>>   sense        : Minimize
>>>   objective    : (x-1)^2+(y-2.5)^2
>>>   #constraints : 3
>>>   #variables   : 2 (Continuous 2)
>>> 
>>>   C 0, name None, -2-(x-(2*y)) <= 0
>>>   C 1, name None, -x <= 0
>>>   C 2, name None, -y <= 0
>>> 
>>>   V 0, name x, Continuous None <= x <= None
>>>   V 1, name y, Continuous None <= y <= None

ラグランジュ関数の生成

まず制約の個数だけ変数を生成します。

# create lagrangian
u = flopt.Variable.array("u", len(prob.constraints), lowBound=0, ini_value=0)

その後、目的関数に係数変数uをかけた制約を足し合わせてラグランジュ関数を作成します。制約はすべて(式) == 0 または (式) <= 0の形で管理されており、.expression属性で取り出すことができます。

const_expressions = [const.expression for const in prob.constraints]
lagrangian = prob.obj + flopt.dot(u, const_expressions)

双対問題の作成

ここで双対問題の目的関数を作成します。双対問題の目的関数も $\min_x L(\theta, x) = f(x) + \sum_i \theta_i g_i(x)$ という最適化問題になっています。そこでCustomExpressionを用いて作成しましょう。

# variables of primary problem
p = prob.getVariables()

# objective function
def dual_obj(u):
    return lagrangian.minimize(optimized_variables=p, timelimit=5)

optimized_variablesによって最適化する変数を制限することができます。ここでは主問題に含まれる変数のみ最適化させたい（すなわち、係数変数uは変化させたくない）ので、それを指定します。そして作成したdual_objを用いて双対問題を定義します。

dual_prob = flopt.Problem(name="Dual problem", sense="Maximize")
dual_prob += flopt.CustomExpression(dual_obj, args=[u])

あとはこれを解くだけです。問題を定義できましたが、Problem::solve()を実行すると、floptの中身では、(1) 主問題の変数(p)を固定してuを最適化 (2) uを固定してpを最適化、が繰り返されることになります。なので双対問題に最適解が存在したとしても、それが得られる保証がないことに注意です。ただ、双対問題であれば最適解が得られなくても主問題の目的関数に対する下限を得ることは可能です。

dual_prob.solve(optimized_variables=u, timelimit=10, msg=True)
print("u", flopt.value(u))
>>> u [0.7999999918945726 0.0 0.0]

解いた結果はこのような感じ。このuを固定してラグランジュ関数の最小化を行います。

lagrangian.minimize(optimized_variables=p, timelimit=10, msg=True)
print("dual obj", lagrangian.value())
>>> dual obj 0.8000000000000002

一連の流れで、主問題の下限が得られました。今回は主問題の目的関数と等しい0.8が得られましたね。（数値誤差の影響で少し大きくなっていますが）

線形計画法の双対問題

線形計画問題の場合はよりタイトな双対問題を作成することができます。

主問題

max  c^Tx
s.t. Ax = b
     x >= 0

について、双対問題を

min  b^Ty
s.t. A^Ty >= c

とすることができます。線形計画問題は非常に良い性質を持っており、もし主問題に最適解が存在すれば、この双対問題にも最適解が存在し、その最適値が一致します（強双対定理）。

線形計画問題probの変数がすべて非負のboundで構成されている場合、双対問題を生成するのは次のコードになります。

import flopt.convert

def dual_lp(prob):
    prob = prob.clone()
    if prob.sense.lower() == "Minimize":
        prob += -prob.obj
    prob = prob.toEq()
    lp = flopt.convert.LpStructure.fromFlopt(prob)
    dual_lp = flopt.convert.LpStructure(c=lp.b, C=lp.C, G=-lp.A.T, h=-lp.c)
    dual_prob = dual_lp.toFlopt(var_name="y")
    dual_prob.sense = "Minimize"
    return dual_prob

Sabrou-mal サブロウ丸

主にプログラミングと数学