数理モデルの定式化#
前節までの説明を踏まえ、本節ではいよいよ数理モデルの定式化方法について述べます。 決定変数やプレースホルダーについては「変数の定義」で触れていますので、本節では目的関数と制約条件の設定方法について説明します。
import jijmodeling as jm
Tip
説明の便宜上、以下では目的関数→制約条件の順に扱いますが、実際のコードでは任意の順番で更新を行うことができます。
目的関数の設定#
Problemオブジェクトの生成時に sense を MAXIMIZE にすると目的関数を最大化する問題、 sense を MINIMIZE にすると最小化する問題として解釈されます。
Problem オブジェクトが作成された初期段階では目的関数は \(0\) として設定され、Problemオブジェクトに対し += 演算子を使って目的関数の項を足していく形で設定します。
Problemオブジェクトが目的関数の項として受け付けるのは、数値型の Expressionオブジェクトのみです。
配列型や辞書型などの式を足そうとすると型エラーとなるので注意してください。
具体的な例として、ナップサック問題の目的関数を設定してみましょう。
@jm.Problem.define(
"Knapsack Problem",
sense=jm.ProblemSense.MAXIMIZE, # 最大化問題として指定
)
def knapsack_problem(problem: jm.DecoratedProblem):
N = problem.Length(description="Number of items")
x = problem.BinaryVar(
shape=(N,), description="$x_i = 1$ if item i is put in the knapsack"
)
v = problem.Float(shape=(N,), description="value of each item")
w = problem.Float(shape=(N,), description="weight of each item")
W = problem.Float(description="maximum weight capacity of the knapsack")
# `+=` 演算子に `Expression` オブジェクトを与えることで目的関数が設定できる
problem += jm.sum(v[i] * x[i] for i in N)
# あるいは、ブロードキャストを用いて次のように書いても「同値」
# problem += jm.sum(v * x)
knapsack_problem
また、JijModeling では、目的関数に項を追加することはできても、全体を書き換えたり削除したりすることはできません。
特に、+= による目的関数の「追加」は新たな項の「追加」として振る舞い、既存の項を別の項で置き換えるものではありません。
次の例を考えます。ここではまず、\(x\)のみを項に持つ目的関数を設定しています。
problem = jm.Problem("Sample")
x = problem.BinaryVar("x")
problem += x
problem
更に、新たな決定変数 \(y\) を定義し、目的関数に \(y\) を追加してみましょう。
y = problem.BinaryVar("y")
problem += y
problem
既存の項が置き換えられたのではなく、\(y\) が加算され \(x + y\) が新たな目的関数となっていることが分かります。 目的関数の項を削除したい場合、目的関数の項の一覧を(Python の)リストなどで持っておき、あとからそれを使って目的関数を設定するなどするとよいでしょう。
制約条件の設定#
制約条件の追加も同様に += 演算子を使って行います。
ただし、制約条件の追加の際には、Problem.Constraint() 関数を使って生成された Constraint オブジェクトを足し合わせる形で追加します。
Problem.Constraint() は必須引数として名前と、==、<=、または >= のいずれかで書かれた制約条件の式を受け取ります。
重要
制約条件の構築に使える比較演算子は ==、<=、>= のみです。
次に示すような > や <、あるいは論理演算などはサポートされていませんので注意してください。
problem.Constraint("BAD1", 1 < x) # ERROR! `>` は使えない!
problem.Constraint("BAD2", (x + y) <= 1 or (y + z) >= 2) # ERROR! 論理演算は使えない!
problem.Constraint("BAD2", (x + y) <= 1 | (y + z) >= 2) # ERROR! 論理演算は使えない!
上で作成したナップサック問題のモデルに制約条件を追加し、モデルを完成させてみましょう。
@knapsack_problem.update
def _(problem: jm.DecoratedProblem):
N = problem.placeholders["N"]
w = problem.placeholders["w"]
W = problem.placeholders["W"]
x = problem.decision_vars["x"]
problem += problem.Constraint("weight", jm.sum(w[i] * x[i] for i in N) <= W)
knapsack_problem
制約条件の追加時には必ず += を呼ぶこと!
制約条件を追加する際には、必ず += 演算子を使って追加してください。単純に Problem.Constraint() を呼び出しただけでは、制約条件はモデルに追加されません。
制約条件の族#
更に、JijModeling では単体の制約条件だけではなく、複数の制約条件をまとめて「制約条件の族」として追加することもできます。 これには複数の方法があります:
domain=や内包表記を使った添え字つき制約条件の定義配列に対する比較式
これらの用法を見るために、ここでは巡回セールスマン問題の二次定式化を考えてみましょう。 都市\(i, j\)の間の距離行列\(d_{i,j}\)、時刻\(t\)に都市\(i\)を訪問することを表すバイナリ変数\(x_{t,i}\)を用いて以下のように表される定式化です:
二種類の制約条件が設定されていますが、それぞれ単一の制約ではなく \(t\) と \(i\) というパラメータを渡る族として定式化されていることに注意しましょう。
添え字つき制約条件#
このように、添え字つきの制約条件を Decorator API で定義するには、Problem.Constraint() メソッドの第二引数をリスト内包表記またはジェネレータ式によって与えればよいです:
@jm.Problem.define("TSP, Decorated", sense=jm.ProblemSense.MINIMIZE)
def tsp_decorated(problem: jm.DecoratedProblem):
C = problem.CategoryLabel(description="Labels of Cities")
N = C.count()
x = problem.BinaryVar(
dict_keys=(N, C), description="$x_{t,i} = 1$ if City $i$ is visited at time $t$"
)
d = problem.Float(dict_keys=(C, C), description="distance between cities")
problem += jm.sum(
d[i, j] * x[t, i] * x[(t + 1) % N, j] for t in N for i in C for j in C
)
# リスト内包表記を使った定義
problem += problem.Constraint(
"one time", [jm.sum(x[t, i] for t in N) == 1 for i in C]
)
# ジェネレータ式を使った定義
problem += problem.Constraint(
"one city", (jm.sum(x[t, i] for i in C) == 1 for t in N)
)
tsp_decorated
Plain API のみで記述する場合は、次のようにパラメータを受け取る lambda 式を第二引数に与え、domain= キーワード引数を合わせて指定します:
tsp_plain = jm.Problem("TSP, Decorated", sense=jm.ProblemSense.MINIMIZE)
C = tsp_plain.CategoryLabel("C", description="Labels of Cities")
N = C.count()
x = tsp_plain.BinaryVar(
"x",
dict_keys=(N, C),
description="$x_{t,i} = 1$ if City $i$ is visited at time $t$",
)
d = tsp_plain.Float("d", dict_keys=(C, C), description="distance between cities")
tsp_plain += jm.sum(
jm.product(N, C, C), lambda t, i, j: d[i, j] * x[t, i] * x[(t + 1) % N, j]
)
# `domain=` キーワード引数を使った定義
tsp_plain += tsp_plain.Constraint(
"one time", lambda i: jm.sum(N, lambda t: x[t, i]) == 1, domain=C
)
tsp_plain += tsp_plain.Constraint(
"one city", lambda t: jm.sum(C, lambda i: x[t, i]) == 1, domain=N
)
tsp_plain
配列同士の比較#
もう一つの方法は、配列や集合の間の比較式を用いて制約条件の族を定義する方法です。 式の構築 で触れたように、比較式にもブロードキャストを用いることができます。 具体的には、制約条件の構築に使える比較式は、両辺が以下の組み合わせのものです:
集合とスカラーの比較
同一シェイプの配列同士の比較
同一キー集合の
TotalDict同士の比較
これを使えば、次のように巡回セールスマン問題の制約条件を定義することができます:
@jm.Problem.define("TSP, Decorated", sense=jm.ProblemSense.MINIMIZE)
def tsp_array_comparison(problem: jm.DecoratedProblem):
N = problem.Natural(description="Number of cities")
x = problem.BinaryVar(
shape=(N, N), description="$x_{t,i} = 1$ if City $i$ is visited at time $t$"
)
d = problem.Float(shape=(N, N), description="distance between cities")
problem += jm.sum(
d[i, j] * x[t, i] * x[(t + 1) % N, j] for t in N for i in N for j in N
)
# 集合とスカラーの比較を使った定義
problem += problem.Constraint("one time", x.sum(axis=0) == 1)
problem += problem.Constraint("one city", x.sum(axis=1) == 1)
tsp_array_comparison
ここで、Expression.sum() や jm.sum() メソッドに axis=i 引数を与えると、Numpy の numpy.sum() 関数と同様に、単純な総和ではなく、その軸に沿った和を計算した配列を返します(複数の軸をリストとして指定することもできます)。
このため、上の例の one-city では x.sum(axis=1) は(\(0\)起点なので)都市を表す\(2\)番目の軸に沿って和を取り、各時刻に訪問される都市の数を表す配列を計算させています。
実際に型を推論させてみると、一次元配列になっているのがわかります。
tsp_array_comparison.infer(tsp_array_comparison.decision_vars["x"].sum(axis=1))
このようにして得られた「時刻毎の都市数」の一元配列をスカラー値\(1\)と比較し、制約条件の族を定義しているのです。 one-timeも同様です。
ここでは配列対スカラーの比較になっていますが、前述の通り同一シェイプの配列同士の比較による制約条件の族の定義も可能です。
同名の制約条件#
複雑な数理モデルを記述していると、場合によっては同名の制約条件を複数回に分けて追加したい場合があります。 JijModeling では以下の条件すべてが満たされている場合に限って、同名の制約条件を複数回に分けて追加することが可能です:
同名の制約条件の添え字の数が一致していること
添え字つきでない制約条件については、添え字の数が \(0\) と見なされます
添え字が重複している場合、制約の定義が description や定義域も含めて厳密に一致していること
(2) の条件については、コンパイル時に与えられたインスタンスデータによって結果がかわってくるため、こうした検査はモデル構築時ではなくインスタンスへのコンパイル時に行われます。これは、添え字の重複の有無は原理的にコンパイル時にしか確定しないためです。次の例を考えてみましょう:
@jm.Problem.define("Possibly Overlapping Constraints")
def problem(problem: jm.DecoratedProblem):
N = problem.Natural(ndim=1)
M = problem.Natural(ndim=1)
n = jm.max(N.max(), M.max()) + 1
x = problem.IntegerVar(shape=n, lower_bound=0, upper_bound=10)
problem += problem.Constraint("constr", [x[i] >= 1 for i in N])
problem += problem.Constraint("constr", [x[i] <= 2 for i in M])
problem
この例では、 \(N\), \(M\) は一次元の添え字の集合であり、データを与えるまで重複しているかどうかはわかりません。たとえば、次のようなインスタンスデータを与えれば、上の問題は問題なくインスタンスへとコンパイルできます:
instance_ok = problem.eval({"N": [0, 2, 4], "M": [1, 3, 5]})
instance_ok.constraints_df
| equality | type | used_ids | name | subscripts | description | |
|---|---|---|---|---|---|---|
| id | ||||||
| 0 | <=0 | Linear | {0} | constr | [0] | <NA> |
| 1 | <=0 | Linear | {2} | constr | [2] | <NA> |
| 2 | <=0 | Linear | {4} | constr | [4] | <NA> |
| 3 | <=0 | Linear | {1} | constr | [1] | <NA> |
| 4 | <=0 | Linear | {3} | constr | [3] | <NA> |
| 5 | <=0 | Linear | {5} | constr | [5] | <NA> |
一方、次のように \(N\) と \(M\) の要素に重複がある場合、コンパイル時エラーとなります:
try:
instance_ng = problem.eval({"N": [0, 2, 4], "M": [1, 2, 5]})
except jm.ModelingError as e:
print(e)
Traceback (most recent last):
while evaluating problem `Problem(name="Possibly Overlapping Constraints", sense=MINIMIZE, objective=0, constraints={constr: [Constraint(name="constr", , lambda i: x[i] >= 1, domain=set(N)),Constraint(name="constr", , lambda i: x[i] <= 2, domain=set(M)),],})',
defined at File "/tmp/ipykernel_715/316836703.py", line 1, col 2-55
File "/tmp/ipykernel_715/316836703.py", line 1, col 2-55:
1 | @jm.Problem.define("Possibly Overlapping Constraints")
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Constraint 'constr' already has subscript 2 with conflicting definion!
old: UserInput(Constraint(name="constr", , lambda i: x[i] >= 1, domain=set(N)))
new: UserInput(Constraint(name="constr", , lambda i: x[i] <= 2, domain=set(M)))
一方、添え字の次元の不一致や、添え字を持たない制約条件の衝突など、直ちに衝突していることがわかるケースは、コンパイル時ではなくモデル構築の時点でエラーとなるようになっています。
try:
@jm.Problem.define("Trivially Conflicting Constraints")
def problem(problem: jm.DecoratedProblem):
x = problem.IntegerVar(lower_bound=0, upper_bound=10)
problem += problem.Constraint("constr", x >= 1)
problem += problem.Constraint("constr", x <= 2)
except jm.ModelingError as e:
print(e)
Traceback (most recent last):
while adding constraint 'constr',
defined at File "/tmp/ipykernel_715/2812696639.py", line 7, col 9-56
File "/tmp/ipykernel_715/2812696639.py", line 7, col 9-56:
7 | problem += problem.Constraint("constr", x <= 2)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Constraint 'constr' has conflicting definition!
existing: Constraint(name="constr", sense=GREATER_THAN_EQUAL, left=x, right=1, shape=Scalar(Integer))
defined at: File "/tmp/ipykernel_715/2812696639.py", line 6, col 20-56
new: Constraint(name="constr", sense=LESS_THAN_EQUAL, left=x, right=2, shape=Scalar(Integer))
defined at: File "/tmp/ipykernel_715/2812696639.py", line 7, col 20-56
とはいえ、制約条件の名前の衝突はコンパイル時にしか確定できないケースがあるため、基本的に同名の制約条件を避けることをお勧めします。