制約最適化システム SCOP （デモ）

Coonstraint Programming Solver SCOP

Open In Colab

SCOP（Solver forCOnstraint Programing：スコープ）は， 大規模な制約最適化問題を高速に解くためのソルバーである．

ここで，制約最適化(constraint optimization)とは， 数理最適化を補完する最適化理論の体系であり， 組合せ最適化問題に特化した求解原理-メタヒューリスティクス(metaheuristics)-を用いるため， 数理最適化ソルバーでは求解が困難な大規模な問題に対しても，効率的に良好な解を探索することができる．

SCOPのトライアルバージョンは， http://logopt.com/scop2/ からダウンロード (Windowsのみ）もしくは

pip install scop

でインストールできる（Pythonのバージョンは3.10以上）。

また，テクニカルドキュメントは，https://scmopt.github.io/moai-manual/scop.html にある．

#Google Colab.では以下でインストールする。
!pip install scop
Requirement already satisfied: scop in /Users/mikiokubo/Library/Caches/pypoetry/virtualenvs/scmopt-HY5JLThu-py3.8/lib/python3.8/site-packages (0.0.1)
WARNING: You are using pip version 20.2.2; however, version 24.3.1 is available.
You should consider upgrading via the '/Users/mikiokubo/Library/Caches/pypoetry/virtualenvs/scmopt-HY5JLThu-py3.8/bin/python -m pip install --upgrade pip' command.

重み付き制約充足問題

ここでは，SCOPで対象とする重み付き制約充足問題について解説する．

一般に制約充足問題(constraint satisfaction problem)は，以下の3つの要素から構成される．

  • 変数(variable): 分からないもの，最適化によって決めるもの． 制約充足問題では，変数は，与えられた集合（以下で述べる「領域」）から1つの要素を選択することによって決められる．

  • 領域(domain): 変数ごとに決められた変数の取り得る値の集合

  • 制約(constraint): 幾つかの変数が同時にとることのできる値に制限を付加するための条件． SCOPでは線形制約（線形式の等式，不等式），2次制約（一般の2次式の等式，不等式）， 相異制約（集合に含まれる変数がすべて異なることを表す制約）が定義できる．

制約充足問題は，制約をできるだけ満たすように， 変数に領域の中の1つの値を割り当てることを目的とした問題である．

SCOPでは，重み付き制約充足問題(weighted constraint satisfaction problem) を対象とする．

ここで「制約の重み」とは，制約の重要度を表す数値であり， SCOPでは正数値もしくは無限大を表す文字列 ’inf’を入力する． ’inf’を入力した場合には，制約は絶対制約(hard constraint)とよばれ， その逸脱量は優先して最小化される． 重みに正数値を入力した場合には，制約は考慮制約(soft constraint)とよばれ， 制約を逸脱した量に重みを乗じたものの和の合計を最小化する．

すべての変数に領域内の値を割り当てたものを解(solution)とよぶ． SCOPでは，単に制約を満たす解を求めるだけでなく， 制約からの逸脱量の重み付き和（ペナルティ）を最小にする解を探索する．

SCOPの基本クラス

SCOPモジュール (scop.py) は，以下のクラスから構成されている．

  • モデルクラス Model

  • 変数クラス Variable

  • 制約クラス Constraint (これは，以下のクラスのスーパークラスである．）

    • 線形制約クラス Linear
    • 2次制約クラス Quadratic
    • 相異制約クラス Alldiff

例題

ここでは，幾つかの簡単な例題を通してSCOPの基本的な使用法を解説する．

以下の例題を動かすためには，最初に以下を追加する必要がある．

from scop import *

もしくは（Marimoのように * でインポートができない場合）

from scop import Model, Variable, Linear, Quadratic, Alldiff
from scop import Model, Variable, Linear, Quadratic, Alldiff

例題 仕事の割当1

あなたは，土木事務所の親方だ．いま，3人の作業員 A,B,C を3つの仕事 \(0,1,2\) に割り当てる必要がある． すべての仕事には1人の作業員を割り当てる必要があるが， 作業員と仕事には相性があり，割り当てにかかる費用（単位は万円）は，以下のようになっているものとする．

仕事 0 1 2
作業員
A 15 20 30
B 7 15 12
C 25 10 13

総費用を最小にするように作業員に仕事を割り振るには，どのようにしたら良いだろうか？

この問題をSCOPを使って解いてみよう．

まず，モデルのインスタンスmodelを生成し，作業員 \(A,B,C\) に割り振られた仕事を表す変数 \(X_A, X_B, X_C\) （プログラムでは A,B,C）を定義する． 数理最適化においては変数は数字で表さなければならないが， 制約最適化では，変数のとれる値の集合で定義する． これを 領域 (domain)とよぶ． 作業員は仕事 \(0,1,2\) のいずれかの仕事をすることができるので，各変数の領域は[0,1,2]となる． 変数の追加は，モデルインスタンスmodelのaddVariableメソッドを用いる．

model = Model()
A = model.addVariable(name="A", domain=[0,1,2])
B = model.addVariable(name="B", domain=[0,1,2])
C = model.addVariable(name="C", domain=[0,1,2])

数理最適化でモデリングをすると \(9=(3 \times 3)\) 個の0-1変数が必要になるが， SCOPだと3個の変数で表現できる．

すべての仕事に1人の作業員を割り当てることを表すには， 相異制約を使う．

  • 相異制約 (Alldiff): リストに含まれる変数（すべて同じ領域をもつと仮定する）がすべて異なる値をとることを表す．

Alldiffクラスのインスタンスalldiffを生成し， それをmodelに追加する． SCOPにおける制約は， すべて逸脱したときのペナルティをweight引数で定義する．weight引数に無限大を表すinfを入れると，絶対制約（ハードな制約）を定義できる．

alldiff = Alldiff("All Diff",[A,B,C],weight="inf")
model.addConstraint(alldiff)

これも数理最適化でモデリングすると，仕事ごとに定義する必要があるので 3本の制約が必要であるが， SCOPだと相異制約1本で表現できる．

SCOPには目的関数という概念がない． すべて制約で表現し， 制約の逸脱ペナルティの合計を最小化する． 割り当て費用は線形制約で記述する．

線形制約は線形不等式（もしくは等式）であり，式として記述する際には，値変数の概念を用いる． 値変数とは変数が領域の値をとったときに \(1\) になる仮想の変数であり， 実際のプログラム内では使わない． 作業員 \(A\) に割り当てられた仕事を表す変数 \(X_A\) に対して，値変数 \(x_{Aj} (j=0,1,2)\) が定義される． \(x_{Aj}\) は， 作業員 \(A\) が仕事 \(j\) に割り当てられたときに \(1\)，それ以外のとき \(0\) を表す変数である．

これを使うと割り当て費用を表す線形制約は，

\[ 15 x_{A0} + 20 x_{A1} + 30 x_{A2} + 7 x_{B0} + 15 x_{B1} + 12 x_{B2} + 25 x_{C0} + 10 x_{C1} + 13 x_{C2} \leq 0 \]

と書ける． この制約の逸脱ペナルティweightを \(1\) に設定すると，制約の逸脱を許す考慮制約（ソフトな制約）となり，逸脱量が割り当て費用になる．

線形制約クラスLinearの右辺定数rhsを \(0\)，制約の方向を<=と設定してインスタンスlinearを生成する． 左辺の各項は，addTermsメソッドを用いて追加する．引数は順に，係数のリスト，変数のリスト，値のリストである．

linear = Linear(name="Objective Function",weight=1,rhs=0,direction="<=")
linear.addTerms([15,20,30],[A,A,A],[0,1,2])
linear.addTerms([7,15,12],[B,B,B],[0,1,2])
linear.addTerms([25,10,13],[C,C,C],[0,1,2])                  
model.addConstraint(linear)

SCOPのインスタンスはprint関数で表示できる． ここでは上で作成したモデルインスタンスmodelを表示しておく． modelのoptimizeメソッドで最適化を実行する．返値は解を表す辞書と逸脱した制約を表す辞書である．

print(model)
sol, violated = model.optimize()
print("solution=", sol)
print("violated constraint=", violated)

プログラム全体を以下に示す．

#単純に記述したプログラム
model = Model()
#変数の宣言
A = model.addVariable(name="A", domain=[0,1,2])
B = model.addVariable(name="B", domain=[0,1,2])
C = model.addVariable(name="C", domain=[0,1,2])

#相異制約
alldiff = Alldiff("All Diff",[A,B,C],weight="inf")
model.addConstraint(alldiff)

#目的関数
linear = Linear(name="Objective Function",weight=1,rhs=0,direction="<=")
linear.addTerms([15,20,30],[A,A,A],[0,1,2])
linear.addTerms([7,15,12],[B,B,B],[0,1,2])
linear.addTerms([25,10,13],[C,C,C],[0,1,2])                  
model.addConstraint(linear)

print(model)

#返値は解を表す辞書と逸脱を表す辞書
sol, violated = model.optimize()
print("solution=", sol)
print("violated constraint=", violated)
Model: 
number of variables = 3  
number of constraints= 2  
variable A:['0', '1', '2'] = None 
variable B:['0', '1', '2'] = None 
variable C:['0', '1', '2'] = None 
All_Diff: weight= inf type=alldiff  A B C ;  :LHS =0  
Objective_Function: weight= 1 type=linear 15(A,0) 20(A,1) 30(A,2) 7(B,0) 15(B,1) 12(B,2) 25(C,0) 10(C,1) 13(C,2) <=0 :LHS =0 

 ================ Now solving the problem ================ 

solution= {'A': '0', 'B': '2', 'C': '1'}
violated constraint= {'Objective_Function': 37}
'''
Example 1 (Assignment Problem):
Three jobs (0,1,2) must be assigned to three workers (A,B,C)
so that each job is assigned to exactly one worker.
The cost matrix is represented by the list of lists
Cost=[[15, 20, 30],
      [7, 15, 12],
      [25,10,13]],
where rows of the matrix are workers, and columns are jobs.
Find the minimum cost assignment of workers to jobs.
'''
#リストと辞書を用いたプログラム
workers=['A','B','C']
Jobs   =[0,1,2]
Cost={ ('A',0):15, ('A',1):20, ('A',2):30,
       ('B',0): 7, ('B',1):15, ('B',2):12,
       ('C',0):25, ('C',1):10, ('C',2):13 }

model=Model()
x={}
for i in workers:
    x[i]=model.addVariable(name=i,domain=Jobs)

xlist=[]
for i in x:
    xlist.append(x[i])

con1=Alldiff('AD',xlist,weight='inf')

con2=Linear('linear_constraint',weight=1,rhs=0,direction='<=')
for i in workers:
    for j in Jobs:
        con2.addTerms(Cost[i,j],x[i],j)

model.addConstraint(con1)
model.addConstraint(con2)

print(model)

model.Params.TimeLimit=1
sol,violated=model.optimize()

if model.Status==0:
    print('solution')
    for x in sol:
        print (x,sol[x])
    print ('violated constraint(s)')
    for v in violated:
        print (v,violated[v])
Model: 
number of variables = 3  
number of constraints= 2  
variable A:['0', '1', '2'] = None 
variable B:['0', '1', '2'] = None 
variable C:['0', '1', '2'] = None 
AD: weight= inf type=alldiff  A C B ;  :LHS =0  
linear_constraint: weight= 1 type=linear 15(A,0) 20(A,1) 30(A,2) 7(B,0) 15(B,1) 12(B,2) 25(C,0) 10(C,1) 13(C,2) <=0 :LHS =0 

 ================ Now solving the problem ================ 

solution
A 0
B 2
C 1
violated constraint(s)
linear_constraint 37

例題 仕事の割当2

あなたは土木事務所の親方だ．今度は，5人の作業員 A,B,C,D,Eを3つの仕事 \(0,1,2\) に割り当てる必要がある． ただし，各仕事にかかる作業員の最低人数が与えられており，それぞれ \(1,2,2\)人必要であり， 割り当ての際の費用（単位は万円）は，以下のようになっているものとする．

仕事 0 1 2
作業員
A 15 20 30
B 7 15 12
C 25 10 13
D 15 18 3
E 5 12 17

さて，誰にどの仕事を割り振れば費用が最小になるだろうか？

例題1では変数をA,B,Cと別々に定義したが，ここではより一般的な記述法を示す．

パラメータは例題1と同じようにリストと辞書で準備する．

  • \(W\)，: 作業員の集合． その要素を \(i\) とする．
  • \(J\): 仕事の集合． その要素を \(j\) とする．
  • \(c_{ij}\): 作業員 \(i\) が仕事 \(j\) に割り当てられたときの費用
  • \(LB_j\): 仕事 \(j\) に必要な人数
model=Model()
workers=['A','B','C','D','E']
Jobs   =[0,1,2]
Cost={ ('A',0):15, ('A',1):20, ('A',2):30,
       ('B',0): 7, ('B',1):15, ('B',2):12,
       ('C',0):25, ('C',1):10, ('C',2):13,
       ('D',0):15, ('D',1):18, ('D',2): 3,
       ('E',0): 5, ('E',1):12, ('E',2):17
       }
LB={0: 1,
    1: 2,
    2: 2
    }

変数は辞書xに保管する．

  • \(X_{i}\): 作業員 \(i\) に割り振られた仕事を表す変数． 領域は仕事の集合 \(J\) であり，そのうち1つの「値」を選択する．
x={}
for i in workers:
    x[i]=model.addVariable(name=i,domain=Jobs)

\(x_{ij}\) は， \(X_i\) が \(j\) に割り当てられたときに \(1\)，それ以外のとき \(0\) を表す変数（値変数）であり，これを使うと人数の下限制約と割り当て費用は，以下の 線形制約として記述できる．

人数下限を表す線形制約（重み \(\infty\)） \[ \sum_{i \in W} x_{ij} \geq LB_j \ \ \ \forall j \in J \]

割り当て費用を表す線形制約（重み \(1\)） \[ \sum_{i \in W, j \in J} c_{ij} x_{ij} \leq 0 \]

LBC={} 
for j in Jobs:
    LBC[j]=Linear(f"LB{j}","inf",LB[j],">=")
    for i in workers:
        LBC[j].addTerms(1,x[i],j)
    model.addConstraint(LBC[j])
    
obj=Linear("obj")
for i in workers:
    for j in [0,1,2]:
        obj.addTerms(Cost[i,j],x[i],j)
model.addConstraint(obj)

modelのパラメータParamsで制限時間TimeLimitを1（秒）に設定して最適化する．

プログラム全体を以下に示す．

"""
Example 2 (Generalized Assignment Problem):
Three jobs (0,1,2) must be assigned to five workers (A,B,C,D,E).
The numbers of workers that must be assigned to jobs 0,1 and 2 are 1,1 and 2, respectively. 
The cost matrix is represented by the list of lists  
Cost=[[15, 20, 30],
[7, 15, 12],
[25,10,13],
[15,18,3],
[5,12,17]]
where rows are workers, and columns are jobs.
Find the minimum cost assignment of workers to jobs.
"""

model=Model()
workers=['A','B','C','D','E']
Jobs   =[0,1,2]
Cost={ ('A',0):15, ('A',1):20, ('A',2):30,
       ('B',0): 7, ('B',1):15, ('B',2):12,
       ('C',0):25, ('C',1):10, ('C',2):13,
       ('D',0):15, ('D',1):18, ('D',2): 3,
       ('E',0): 5, ('E',1):12, ('E',2):17
       }
LB={0: 1,
    1: 2,
    2: 2
    }
x={}
for i in workers:
    x[i]=model.addVariable(name=i,domain=Jobs)
LBC={} 
for j in Jobs:
    LBC[j]=Linear(f"LB{j}","inf",LB[j],">=")
    for i in workers:
        LBC[j].addTerms(1,x[i],j)
    model.addConstraint(LBC[j])
obj=Linear("obj")
for i in workers:
    for j in [0,1,2]:
        obj.addTerms(Cost[i,j],x[i],j)
model.addConstraint(obj)

model.Params.TimeLimit=1
sol,violated=model.optimize()

print('solution')
for x in sol:
    print (x,sol[x])
print ('violated constraint(s)')
for v in violated:
    print (v,violated[v])

 ================ Now solving the problem ================ 

solution
A 1
B 2
C 1
D 2
E 0
violated constraint(s)
obj 50

例題 仕事の割当3

上の例題と同じ状況で，仕事を割り振ろうとしたところ，作業員 A と C は仲が悪く， 一緒に仕事をさせると喧嘩を始めることが判明した． 作業員 A と C を同じ仕事に割り振らないようにするには，どうしたら良いだろうか？

この問題は，追加された作業員 A と C を同じ仕事に割り当てることを禁止する制約を記述するだけで解決できる． ここでは，2次制約（重みは \(100\)）として記述する．

\[ x_{A0} x_{C0} + x_{A1} x_{C1} + x_{A2} x_{C2} = 0 \]

作業員AとCが同じ仕事に割り当てられると左辺は \(1\)になり，制約を逸脱する．

線形制約クラスと同様に2次制約クラスQuadraticからインスタンスconfを生成する． 左辺の項を追加するには，addTermsメソッドを用いる． 引数は，最初の変数の係数，変数，値の次に2番目の変数の係数，変数，値を入れる．

conf=Quadratic("conflict",100,0,"=")
for j in Jobs:
    conf.addTerms(1,x["A"],j,x["C"],j)
model.addConstraint(conf)

数理最適化ソルバーは非凸の2次を含む制約や目的関数が苦手であるが，SCOPは通常の制約と同じように解くことができる．

"""
Example 3 (Variation of Generalized Assignment Problem):
Three jobs (0,1,2) must be assigned to five workers (A,B,C,D,E).
The minimum numbers of workers that must be assigned to jobs 0,1 and 2 are 1,2 and 2, respectively.
This lower bound is represented by a dictionary:
LB={0: 1,
    1: 2,
    2: 2
    }
where keys are jobs and values are lower bounds.
The cost matrix is represented by a dictionary:
Cost={ ("A",0):15, ("A",1):20, ("A",2):30,
       ("B",0): 7, ("B",1):15, ("B",2):12,
       ("C",0):25, ("C",1):10, ("C",2):13,
       ("D",0):15, ("D",1):18, ("D",2): 3,
       ("E",0): 5, ("E",1):12, ("E",2):17
       }
where keys are tuples of workers and jobs, and values are costs.
We add an additional condition: worker A cannot do the job with worker C.
Find the minimum cost assignment of workers to jobs.
"""

model=Model()
workers=["A","B","C","D","E"]
Jobs   =[0,1,2]
Cost={ ("A",0):15, ("A",1):20, ("A",2):30,
       ("B",0): 7, ("B",1):15, ("B",2):12,
       ("C",0):25, ("C",1):10, ("C",2):13,
       ("D",0):15, ("D",1):18, ("D",2): 3,
       ("E",0): 5, ("E",1):12, ("E",2):17
       }
LB={0: 1,
    1: 2,
    2: 2
    }
x={}
for i in workers:
    x[i]= model.addVariable(i,Jobs)
LBC={}
for j in Jobs:
    LBC[j]=Linear(f"LB{j}","inf",LB[j],">=")
    for i in workers:
        LBC[j].addTerms(1,x[i],j)
    model.addConstraint(LBC[j])
obj=Linear("obj",1,0,"<=")
for i in workers:
    for j in Jobs:
        obj.addTerms(Cost[i,j],x[i],j)
model.addConstraint(obj)
conf=Quadratic("conflict",100,0,"=")
for j in Jobs:
    conf.addTerms(1,x["A"],j,x["C"],j)
model.addConstraint(conf)
model.Params.TimeLimit=1
sol,violated= model.optimize()
print ("solution")
for x in sol:
    print (x,sol[x])
print ("violated constraint(s)")
for v in violated:
    print (v,violated[v])

 ================ Now solving the problem ================ 

solution
A 0
B 2
C 1
D 2
E 1
violated constraint(s)
obj 52

例題 魔方陣

魔方陣とは， \(n \times n\) のマス目に \(1\) から \(n^2\) までの数字を1つずつ入れて，どの横行，縦列，対角線のマス目の数字の和も同じになるようにしたものである.

\(n=3\) の問題を以下の手順で問題を解く．

  1. 各マス目 \((i,j), i=0,1,2, j=0,1,2\) に対して変数 \(x[i,j]\) を準備して，その領域を \(1\) から \(9\) までの数とする．

  2. 各マス目には異なる数字を入れる必要があるので，すべての変数のリストを入れた相異制約 (Alldiff) を追加する． この制約は絶対制約とする．

  3. さらに，各行（\(i=0,1,2\))と各列(\(j=0,1,2\))に対して，その和がちょうど \(15 = (1+2+\cdots+9)/3\) になるという制約を追加する． これらの制約は考慮制約とし，逸脱の重みは \(1\) とする．

  4. 最適化を行い，解を表示する．

  • 行の集合を \(I\)， その要素を \(i\) とする．

  • 列の集合を \(J\)，その要素を \(j\) とする．

  • \(X_{ij}\): マス目 \(i,j\) に割り当てられた数字を表す変数； 領域は \([1,2,3,4,5,6,7,8,9]\) であり，そのうち1つの「値」を選択する．

  • \(x_{ijk}\): \(X_{ij}\) が \(k\) に割り当てられたときに \(1\)，それ以外のとき \(0\) を表す変数（値変数）

相異制約（重み \(\infty\)）； すべてのマス目の数字が異なることを表す．

Alldiff( [ X_{ij} for i in I for j in J  ] )

線形制約（重み \(1\)）；行ごとの和が \(15\) であることを表す． \[ \sum_{j \in J} \sum_{k} k x_{ijk} = 15 \ \ \ \forall i \in I \]

線形制約（重み \(1\)）；列ごとの和が \(15\) であることを表す． \[ \sum_{i \in I} \sum_{k} k x_{ijk} = 15 \ \ \ \forall j \in J \]

線形制約（重み \(1\)）；対角線ごとの和が \(15\) であることを表す． \[ \sum_{j \in J} \sum_{k} k x_{jjk} = 15 \]

\[ \sum_{j \in J} \sum_{k} k x_{j,2-j,k} = 15 \]

以下に一般の \(n\) でも解けるプログラムを示す． ただしトライアル版だと \(n=3\) までしか解くことができない．

n = 3
nn = n*n
model = Model()
x = {}
dom = [i+1 for i in range(nn)]
sum_ = sum(dom)//n 
for i in range(n):
    for j in range(n):
        x[i,j] = model.addVariable(name=f"x[{i},{j}]", domain=dom)
alldiff = Alldiff(f'AD',[ x[i,j] for i in range(n) for j in range(n) ], weight='inf')
model.addConstraint( alldiff )
col_constr = {}
for j in range(n):
    col_constr[j] =  Linear(f'col_constraint{j}',weight=1,rhs=sum_,direction='=')
    for i in range(n):
        for k in range(1,nn+1):
            col_constr[j].addTerms(k,x[i,j],k) 
    model.addConstraint(col_constr[j])
row_constr = {}
for i in range(n):
    row_constr[i] =  Linear(f'row_constraint{i}',weight=1,rhs=sum_,direction='=')
    for j in range(n):
        for k in range(1,nn+1):
            row_constr[i].addTerms(k,x[i,j],k) 
    model.addConstraint(row_constr[i])
diagonal_constr = {}
diagonal_constr[0] =  Linear(f'diagonal_constraint{0}',weight=1,rhs=sum_,direction='=')
for j in range(n):
    for k in range(1,nn+1):
        diagonal_constr[0].addTerms(k,x[j,j],k) 
model.addConstraint(diagonal_constr[0])
diagonal_constr[1] =  Linear(f'diagonal_constraint{1}',weight=1,rhs=sum_,direction='=')
for j in range(n):
    for k in range(1,nn+1):
        diagonal_constr[1].addTerms(k,x[j,n-1-j],k) 
model.addConstraint(diagonal_constr[1])
model.Params.TimeLimit=100
model.Params.RandomSeed=1
#model.Params.OutputFlag=True
sol,violated = model.optimize()
print("逸脱制約=", violated)
import numpy as np
solution = np.zeros( (n,n), int )
for i in range(n):
    for j in range(n):
        solution[i,j] = int(x[i,j].value)
print(solution)

 ================ Now solving the problem ================ 

逸脱制約= {}
[[2 9 4]
 [7 5 3]
 [6 1 8]]

練習問題

問題 多制約ナップサック

あなたは，ぬいぐるみ専門の泥棒だ． ある晩，あなたは高級ぬいぐるみ店にこっそり忍び込んで，盗む物を選んでいる． 狙いはもちろん，マニアの間で高額で取り引きされているクマさん人形だ． クマさん人形は，現在 \(4\)体販売されていて， それらの値段と重さと容積は，以下のリストで与えられている．

v=[16,19,23,28]                     #価値
a=[[2,3,4,5],[3000,3500,5100,7200]] #重さと容積

あなたは，転売価格の合計が最大になるようにクマさん人形を選んで逃げようと思っているが， あなたが逃走用に愛用しているナップサックはとても古く， \(7\)kgより重い荷物を入れると，底がぬけてしまうし，\(10000 {cm}^3\)（\(10\)\(\ell\)）を超えた荷物を入れると破けてしまう．

さて，どのクマさん人形をもって逃げれば良いだろうか？

#
model=Model()

v=[16,19,23,28]
a=[[2,3,4,5],[3000,3500,5100,7200]]
b=[7,10000]
n=len(v)
m=len(b)
items=["item{0}".format(j) for j in range(n)]
varlist=model.addVariables(items,[0,1])
for i in range(m):
    con1=Linear("mkp_{0}".format(i),"inf",b[i])
    for j in range(n):
        con1.addTerms(a[i][j],varlist[j],1)
    model.addConstraint(con1)

con2=Linear("obj",1,sum(v),">=")
for j in range(n):
    con2.addTerms(v[j],varlist[j],1)
model.addConstraint(con2)

model.Params.TimeLimit=1
sol,violated=model.optimize()

print (model)

if model.Status==0:
    print ("solution")
    for x in sol:
        print (x,sol[x])
    print ("violated constraint(s)")
    for v in violated:
        print (v,violated[v])

 ================ Now solving the problem ================ 

Model: 
number of variables = 4  
number of constraints= 3  
variable item0:['0', '1'] = 0 
variable item1:['0', '1'] = 1 
variable item2:['0', '1'] = 1 
variable item3:['0', '1'] = 0 
mkp_0: weight= inf type=linear 2(item0,1) 3(item1,1) 4(item2,1) 5(item3,1) <=7 :LHS =7  
mkp_1: weight= inf type=linear 3000(item0,1) 3500(item1,1) 5100(item2,1) 7200(item3,1) <=10000 :LHS =8600  
obj: weight= 1 type=linear 16(item0,1) 19(item1,1) 23(item2,1) 28(item3,1) >=86 :LHS =42 
solution
item0 0
item1 1
item2 1
item3 0
violated constraint(s)
obj 44

問題 最大安定集合

あなたは \(6\)人のお友達から何人か選んで一緒にピクニックに行こうと思っている． しかし，グラフ上で隣接している（線で結ばれている）人同士はとても仲が悪く，彼らが一緒にピクニックに 行くとせっかくの楽しいピクニックが台無しになってしまう． なるべくたくさんの仲間でピクニックに行くには誰を誘えばいいんだろう？

ただし，グラフの隣接点の情報は以下のリストで与えられているものとする．

adj=[[2],[3],[0,3,4,5],[1,2,5],[2],[2,3]]
#
m=Model()

nodes=["n{0}".format(i) for i in range(6)]
adj=[[2],[3],[0,3,4,5],[1,2,5],[2],[2,3]]
n=len(nodes)

varlist=m.addVariables(nodes,[0,1])

for i in range(n):
    for j in adj[i]:
        if i<j:
            con1=Linear("constraint{0}_{1}".format(i,j),"inf",1)
            con1.addTerms(1,varlist[i],1)
            con1.addTerms(1,varlist[j],1)
            m.addConstraint(con1)

obj=Linear("obj",1,n,">=")
for i in range(n):
    obj.addTerms(1,varlist[i],1)
m.addConstraint(obj)

m.Params.TimeLimit=1
sol,violated=m.optimize()

print (m)

if m.Status==0:
    print ("solution")
    for x in sol:
        print (x,sol[x])
    print ("violated constraint(s)")
    for v in violated:
        print (v,violated[v])

 ================ Now solving the problem ================ 

Model: 
number of variables = 6  
number of constraints= 7  
variable n0:['0', '1'] = 1 
variable n1:['0', '1'] = 1 
variable n2:['0', '1'] = 0 
variable n3:['0', '1'] = 0 
variable n4:['0', '1'] = 1 
variable n5:['0', '1'] = 1 
constraint0_2: weight= inf type=linear 1(n0,1) 1(n2,1) <=1 :LHS =1  
constraint1_3: weight= inf type=linear 1(n1,1) 1(n3,1) <=1 :LHS =1  
constraint2_3: weight= inf type=linear 1(n2,1) 1(n3,1) <=1 :LHS =0  
constraint2_4: weight= inf type=linear 1(n2,1) 1(n4,1) <=1 :LHS =1  
constraint2_5: weight= inf type=linear 1(n2,1) 1(n5,1) <=1 :LHS =1  
constraint3_5: weight= inf type=linear 1(n3,1) 1(n5,1) <=1 :LHS =1  
obj: weight= 1 type=linear 1(n0,1) 1(n1,1) 1(n2,1) 1(n3,1) 1(n4,1) 1(n5,1) >=6 :LHS =4 
solution
n0 1
n1 1
n2 0
n3 0
n4 1
n5 1
violated constraint(s)
obj 2

問題 グラフ彩色

今度は，同じお友達のクラス分けで悩んでいる． お友達同士で仲が悪い組は，グラフ上で隣接している． 仲が悪いお友達を同じクラスに入れると喧嘩を始めてしまう． なるべく少ないクラスに分けるには，どのようにすればいいんだろう？

ただし，グラフの隣接点の情報は以下のリストで与えられているものとする．

adj=[[2],[3],[0,3,4,5],[1,2,5],[2],[2,3]]
#
m=Model()

K=3
nodes=["n{0}".format(i) for i in range(6)]
adj=[[2],[3],[0,3,4,5],[1,2,5],[2],[2,3]]
n=len(nodes)

varlist=m.addVariables(nodes,range(K))

for i in range(n):
    for j in adj[i]:
        if i<j:
            con1=Alldiff("alldiff_{0}_{1}".format(i,j),[varlist[i],varlist[j]],"inf")      
            m.addConstraint(con1)
            
m.Params.TimeLimit=1
sol,violated=m.optimize()

print (m)
if m.Status==0:
    print ("solution")
    for x in sol:
        print (x,sol[x])
    print ("violated constraint(s)")
    for v in violated:
        print (v,violated[v])

 ================ Now solving the problem ================ 

Model: 
number of variables = 6  
number of constraints= 6  
variable n0:['0', '1', '2'] = 0 
variable n1:['0', '1', '2'] = 0 
variable n2:['0', '1', '2'] = 1 
variable n3:['0', '1', '2'] = 2 
variable n4:['0', '1', '2'] = 2 
variable n5:['0', '1', '2'] = 0 
alldiff_0_2: weight= inf type=alldiff  n0 n2 ;  :LHS =0  
alldiff_1_3: weight= inf type=alldiff  n1 n3 ;  :LHS =0  
alldiff_2_3: weight= inf type=alldiff  n3 n2 ;  :LHS =0  
alldiff_2_4: weight= inf type=alldiff  n2 n4 ;  :LHS =0  
alldiff_2_5: weight= inf type=alldiff  n5 n2 ;  :LHS =0  
alldiff_3_5: weight= inf type=alldiff  n5 n3 ;  :LHS =0 
solution
n0 0
n1 0
n2 1
n3 2
n4 2
n5 0
violated constraint(s)

問題 グラフ分割

今度は，同じ\(6\)人のお友達を2つのチームに分けてミニサッカーをしようとしている． もちろん，公平を期すために，同じ人数になるように3人ずつに分ける． ただし，仲が悪いお友達が同じチームになることは極力避けたいと考えている． さて，どのようにチーム分けをしたら良いだろうか？

ただし，中の悪い同士を表すグラフの隣接点の情報は以下のリストで与えられているものとする．

adj=[[1,4],[0,2,4],[1],[4,5],[0,1,3,5],[3,4]]
#
nodes=[f"n{i}" for i in range(6)]
adj=[[1,4],[0,2,4],[1],[4,5],[0,1,3,5],[3,4]]
n=len(nodes)

m = Model()

varlist=m.addVariables(nodes,[0,1])

con1=Linear("constraint","inf",n//2,"=")
for i in range(len(nodes)):
    con1.addTerms(1,varlist[i],1)
m.addConstraint(con1)

##con2={}
##for i in range(n):
##    for j in adj[i]:
##        con2[i,j]= Quadratic( "obj_%s_%s"%(i,j) )
##        con2[i,j].addTerms(1,varlist[i],1,varlist[j],0)
##        con2[i,j].addTerms(1,varlist[i],0,varlist[j],1)
##        m.addConstraint(con2[i,j])

con2=Quadratic( "obj")
for i in range(n):
    for j in adj[i]:
        con2.addTerms(1,varlist[i],1,varlist[j],0)
        con2.addTerms(1,varlist[i],0,varlist[j],1)
m.addConstraint(con2)

print (m)

m.Params.TimeLimit=1
sol,violated=m.optimize()

if m.Status==0:
    print ("solution")
    for x in sol:
        print (x,sol[x])
    print ("violated constraint(s)")
    for v in violated:
        print (v,violated[v])
Model: 
number of variables = 6  
number of constraints= 2  
variable n0:['0', '1'] = None 
variable n1:['0', '1'] = None 
variable n2:['0', '1'] = None 
variable n3:['0', '1'] = None 
variable n4:['0', '1'] = None 
variable n5:['0', '1'] = None 
constraint: weight= inf type=linear 1(n0,1) 1(n1,1) 1(n2,1) 1(n3,1) 1(n4,1) 1(n5,1) =3 :LHS =0  
obj: weight=1 type=quadratic 1(n0,1)(n1,0) 1(n0,0)(n1,1) 1(n0,1)(n4,0) 1(n0,0)(n4,1) 1(n1,1)(n0,0) 1(n1,0)(n0,1) 1(n1,1)(n2,0) 1(n1,0)(n2,1) 1(n1,1)(n4,0) 1(n1,0)(n4,1) 1(n2,1)(n1,0) 1(n2,0)(n1,1) 1(n3,1)(n4,0) 1(n3,0)(n4,1) 1(n3,1)(n5,0) 1(n3,0)(n5,1) 1(n4,1)(n0,0) 1(n4,0)(n0,1) 1(n4,1)(n1,0) 1(n4,0)(n1,1) 1(n4,1)(n3,0) 1(n4,0)(n3,1) 1(n4,1)(n5,0) 1(n4,0)(n5,1) 1(n5,1)(n3,0) 1(n5,0)(n3,1) 1(n5,1)(n4,0) 1(n5,0)(n4,1) <=0 :LHS =0 

 ================ Now solving the problem ================ 

solution
n0 0
n1 0
n2 0
n3 1
n4 1
n5 1
violated constraint(s)
obj 4

問題 巡回セールスマン

あなたは休暇を利用してヨーロッパめぐりをしようと考えている． 現在スイスのチューリッヒに宿を構えているあなたの目的は， スペインのマドリッドで闘牛を見ること， イギリスのロンドンでビックベンを見物すること， イタリアのローマでコロシアムを見ること， ドイツのベルリンで本場のビールを飲むことである．

あなたはレンタルヘリコプターを借りてまわることにしたが， 移動距離に比例した高額なレンタル料を支払わなければならない． したがって， あなたはチューリッヒ (T) を出発した後， なるべく短い距離で他の \(4\)つの都市 マドリッド(M)，ロンドン(L)，ローマ(R)，ベルリン(B) を経由し， 再びチューリッヒに帰って来ようと考えた． 都市の間の移動距離を測ってみたところ，以下のようになっていることがわかった．

cities=["T","L","M","R","B"]
d=[[0,476,774,434,408],
   [476,0,784,894,569],
   [774,784,0,852,1154],
   [434,894,852,0,569],
   [408,569,1154,569,0]]

さて，どのような順序で旅行すれば，移動距離が最小になるだろうか?

#
m=Model()

cities=["T","L","M","R","B"]
d=[[0,476,774,434,408],[476,0,784,894,569],[774,784,0,852,1154],[434,894,852,0,569],[408,569,1154,569,0]]
n=len(cities)

varlist=m.addVariables(cities,range(n))

con1=Alldiff("AD",varlist,"inf")
m.addConstraint(con1)

obj=Quadratic("obj")
for i in range(n):
    for j in range(n):
        if i!=j:
            for k in range(n):
                if k ==n-1:
                    ell=0
                else:
                    ell=k+1
                obj.addTerms(d[i][j],varlist[i],k,varlist[j],ell)
m.addConstraint(obj)

m.Params.TimeLimit=1
sol,violated=m.optimize()

print (m)

if m.Status==0:
    print ("solution")
    for x in sol:
        print (x,sol[x])
    print ("violated constraint(s)")
    for v in violated:
        print (v,violated[v])

 ================ Now solving the problem ================ 

Model: 
number of variables = 5  
number of constraints= 2  
variable T:['0', '1', '2', '3', '4'] = 1 
variable L:['0', '1', '2', '3', '4'] = 3 
variable M:['0', '1', '2', '3', '4'] = 4 
variable R:['0', '1', '2', '3', '4'] = 0 
variable B:['0', '1', '2', '3', '4'] = 2 
AD: weight= inf type=alldiff  B L R M T ;  :LHS =0  
obj: weight=1 type=quadratic 476(T,0)(L,1) 476(T,1)(L,2) 476(T,2)(L,3) 476(T,3)(L,4) 476(T,4)(L,0) 774(T,0)(M,1) 774(T,1)(M,2) 774(T,2)(M,3) 774(T,3)(M,4) 774(T,4)(M,0) 434(T,0)(R,1) 434(T,1)(R,2) 434(T,2)(R,3) 434(T,3)(R,4) 434(T,4)(R,0) 408(T,0)(B,1) 408(T,1)(B,2) 408(T,2)(B,3) 408(T,3)(B,4) 408(T,4)(B,0) 476(L,0)(T,1) 476(L,1)(T,2) 476(L,2)(T,3) 476(L,3)(T,4) 476(L,4)(T,0) 784(L,0)(M,1) 784(L,1)(M,2) 784(L,2)(M,3) 784(L,3)(M,4) 784(L,4)(M,0) 894(L,0)(R,1) 894(L,1)(R,2) 894(L,2)(R,3) 894(L,3)(R,4) 894(L,4)(R,0) 569(L,0)(B,1) 569(L,1)(B,2) 569(L,2)(B,3) 569(L,3)(B,4) 569(L,4)(B,0) 774(M,0)(T,1) 774(M,1)(T,2) 774(M,2)(T,3) 774(M,3)(T,4) 774(M,4)(T,0) 784(M,0)(L,1) 784(M,1)(L,2) 784(M,2)(L,3) 784(M,3)(L,4) 784(M,4)(L,0) 852(M,0)(R,1) 852(M,1)(R,2) 852(M,2)(R,3) 852(M,3)(R,4) 852(M,4)(R,0) 1154(M,0)(B,1) 1154(M,1)(B,2) 1154(M,2)(B,3) 1154(M,3)(B,4) 1154(M,4)(B,0) 434(R,0)(T,1) 434(R,1)(T,2) 434(R,2)(T,3) 434(R,3)(T,4) 434(R,4)(T,0) 894(R,0)(L,1) 894(R,1)(L,2) 894(R,2)(L,3) 894(R,3)(L,4) 894(R,4)(L,0) 852(R,0)(M,1) 852(R,1)(M,2) 852(R,2)(M,3) 852(R,3)(M,4) 852(R,4)(M,0) 569(R,0)(B,1) 569(R,1)(B,2) 569(R,2)(B,3) 569(R,3)(B,4) 569(R,4)(B,0) 408(B,0)(T,1) 408(B,1)(T,2) 408(B,2)(T,3) 408(B,3)(T,4) 408(B,4)(T,0) 569(B,0)(L,1) 569(B,1)(L,2) 569(B,2)(L,3) 569(B,3)(L,4) 569(B,4)(L,0) 1154(B,0)(M,1) 1154(B,1)(M,2) 1154(B,2)(M,3) 1154(B,3)(M,4) 1154(B,4)(M,0) 569(B,0)(R,1) 569(B,1)(R,2) 569(B,2)(R,3) 569(B,3)(R,4) 569(B,4)(R,0) <=0 :LHS =3047 
solution
T 1
L 3
M 4
R 0
B 2
violated constraint(s)
obj 3047

問題 ビンパッキング

あなたは，大企業の箱詰め担当部長だ．あなたの仕事は，色々な大きさのものを，決められた大きさの箱に「上手に」詰めることである． この際，使う箱の数をなるべく少なくすることが，あなたの目標だ． （なぜって，あなたの会社が利用している宅配業者では，運賃は箱の数に比例して決められるから．） 1つの箱に詰められる荷物の上限は \(7\)kgと決まっており，荷物の重さはのリストは [6,5,4,3,1,2] である． しかも，あなたの会社で扱っている荷物は，どれも重たいものばかりなので，容積は気にする必要はない （すなわち箱の容量は十分と仮定する）． さて，どのように詰めて運んだら良いだろうか？

#
bpp=Model()

Items=[6,5,4,3,1,2]
B=7
num_bins=3
n=len(Items)

x={}
for i in range(n):
     x[i] = bpp.addVariable("x_{0}".format(i),range(num_bins))
Bin={}
for j in range(num_bins):
     Bin[j]=Linear("Bin_{0}".format(j),weight=1,rhs=B,direction="<=")
     for i in range(n):
          Bin[j].addTerms(Items[i],x[i],j)
     bpp.addConstraint(Bin[j])

sol,violated=bpp.optimize()

if bpp.Status==0:
    print ("solution=")
    for i in sol:
        print (i,sol[i])

    print ("violated constraints=",violated)

 ================ Now solving the problem ================ 

solution=
x_0 2
x_1 0
x_2 1
x_3 1
x_4 2
x_5 0
violated constraints= {}

問題 最適化版の\(8\)-クイーン

\(8 \times 8\) のチェス盤に \(8\)個のクイーンを置くことを考える． チェスのクイーンとは，将棋の飛車と角の両方の動きができる最強の駒である． \(i\)行 \(j\)列に置いたときの費用を \(i \times j\) と定義したとき， クイーンがお互いに取り合わないように置く配置の中で，費用の合計が最小になるような配置を求めよ．

#
m=Model()

n=8
varlist=[]
for i in range(n):
    varlist.append("x{0}".format(i))

var=m.addVariables(varlist,range(n))

con1=Alldiff("AD",var,"inf")
m.addConstraint(con1)

for k in range(2,2*n-1):
    con2=Linear("rightdown_{0}".format(k),"inf",1,"<=")
    for i in range(n):
        j=k-n+i
        if j>=0 and j<=n-1:
            con2.addTerms(1,var[i],j)
    m.addConstraint(con2)
        
for k in range(2,2*n-1):
    con3=Linear("leftdown_{0}".format(k),"inf",1,"<=")
    for i in range(n):
        j=k-i-1
        if j>=0 and j<=n-1:
            con3.addTerms(1,var[i],j)
    m.addConstraint(con3)

obj=Linear("obj",1,0,"<=")
for i in range(n):
    for j in range(n):
        obj.addTerms((i+1)*(j+1),var[i],j)
m.addConstraint(obj)

m.Params.TimeLimit=1
sol,violated=m.optimize()

#print (m)

if m.Status==0:
    print ("solution")
    for x in sol:
        print (x,sol[x])
    print ("violated constraint(s)")
    for v in violated:
        print (v,violated[v])

 ================ Now solving the problem ================ 

solution
x0 6
x1 3
x2 1
x3 7
x4 5
x5 0
x6 2
x7 4
violated constraint(s)
obj 150

問題 2次割当

いま，3人のお友達が3箇所の家に住もうとしている． 3人は毎週何回か重要な打ち合わせをする必要があり，打ち合わせの頻度は，リストのリスト．

f = [[0,5,1],[5,0,2],[1,2,0]]

として与えられている．

また，家の間の移動距離もリストのリスト

d = [[0,2,3],[2,0,1],[3,1,0]]

として与えられているものとする．

3人は打ち合わせのときに移動する距離を最小に するような場所に住むことを希望している．さて，誰をどの家に割り当てたらよいのだろうか？

#
m=Model()

n=3
d=[[0,2,3],[2,0,1],[3,1,0]]
f=[[0,5,1],[5,0,2],[1,2,0]]

nodes=["n{0}".format(i) for i in range(n)]

varlist=m.addVariables(nodes,range(n))

con1=Alldiff("AD",varlist,"inf")
m.addConstraint(con1)

obj=Quadratic("obj")
for i in range(n-1):
    for j in range(i+1,n):
        for k in range(n):
            for ell in range(n):
                if k !=ell:
                    obj.addTerms(f[i][j]*d[k][ell],varlist[i],k,varlist[j],ell)
m.addConstraint(obj)

m.Params.TimeLimit=1
sol,violated=m.optimize()

if m.Status==0:
    print ("solution")
    for x in sol:
        print (x,sol[x])
    print ("violated constraint(s)")
    for v in violated:
        print (v,violated[v])

 ================ Now solving the problem ================ 

solution
n0 2
n1 1
n2 0
violated constraint(s)
obj 12

問題 車の投入順決定

コンベア上に一直線に並んだ車の生産ラインを考える． このラインは，幾つかの作業場から構成され，それぞれの作業場では異なる作業が行われる． いま，4種類の車を同じ生産ラインで製造しており，それぞれをモデル \(A,B,C,D\) とする． 本日の製造目標は，それぞれ \(30,30,20,40\)台である．

最初の作業場では，サンルーフの取り付けを行っており，これはモデル \(B,C\) だけに必要な作業である． 次の作業場では，カーナビの取り付けが行われており，これはモデル \(A,C\) だけに必要な作業である． それぞれの作業場は長さをもち， サンルーフ取り付けは車 \(5\)台分，カーナビ取り付けは車 \(3\)台分の長さをもつ． また，作業場には作業員が割り当てられており，サンルーフ取り付けは \(3\)人，カーナビ取り付けは \(2\)人の 作業員が配置されており，作業場の長さを超えない範囲で別々に作業を行う．

作業場の範囲で作業が可能な車の投入順序を求めよ．

ヒント： 投入順序をうまく決めないと，作業場の範囲内で作業を完了することができない． たとえば，\(C,A,A,B,C\) の順で投入すると， サンルーフ取り付けでは，3人の作業員がそれぞれモデル \(C,B,C\) に対する作業を行うので 間に合うが，カーナビ取り付けでは， 2人の作業員では \(C,A,A\) の3台の車の作業を終えることができない．

これは，作業場の容量制約とよばれ，サンルーフ取り付けの作業場では， すべての連続する \(5\)台の車の中に，モデル \(B,C\) が高々 \(3\)つ， カーナビ取り付けの作業場では， すべての連続する \(3\)台の車の中に，モデル \(A,C\) が高々 \(2\)つ入っているという制約を課すことに相当する

#
m=Model()
Type=["A","B","C","D","E","F"] #car types
Number={"A":1,"B":1,"C":2,"D":2,"E":2,"F":2}   #number of cars needed 
n=sum(Number[i] for i in Number) #planning horizon
#1st line produces car type B and C that has a workplace with length 5 and 3 workers
#2nd line produces car type A anc C that has a workplace with length 3 and 2 workers 
Option=[["A","E","F"],
    ["C","D","F"],
    ["A","E"],
    ["A","B","D"],
    ["C"]] 
Length=[2,3,3,5,5] 
Capacity=[1,2,1,2,1]

X={}
for i in range(n):
    X[i]=m.addVariable("seq[{0}]".format(i),Type)

#production volume constraints
for i in Type:
    L1=Linear("req[{0}]".format(i),direction="=",rhs=Number[i])
    for j in range(n):
        L1.addTerms(1,X[j],i)
    m.addConstraint(L1)
    
for i in range(len(Length)):
    for k in range(n-Length[i]+1):
        L2=Linear("ub[{0}_{1}]".format(i,k),direction="<=",rhs=Capacity[i])
        for t in range(k,k+Length[i]):
            for j in range(len(Option[i])):
                L2.addTerms(1,X[t],Option[i][j])
        m.addConstraint(L2)

m.Params.TimeLimit=1
m.Params.OutputFlag=False
sol,violated=m.optimize()

if m.Status==0:
    print ("solution")
    for x in sol:
        print (x,sol[x])
    print ("violated constraint(s)")
    for v in violated:
        print (v,violated[v])

 ================ Now solving the problem ================ 

solution
seq[0] A
seq[1] C
seq[2] F
seq[3] B
seq[4] F
seq[5] D
seq[6] E
seq[7] C
seq[8] D
seq[9] E
violated constraint(s)

問題 段取り費用付き生産計画

１つの生産ラインでa,bの２種類の製品を生産している．各期に生産できる製品は1つであり，生産はバッチで行われるため生産量は決まっている（辞書S）． 5期の需要量（辞書D）を満たすように，生産計画（どの期にどの製品を生産するか）を作りたいのだが，製品の切り替えには段取り費用（辞書F）がかかる． ただし，生産しないことを表すダミーの製品０があるものと仮定し，直前の期では何も生産していなかったものと仮定する． 生産すると生産量だけ在庫が増え，毎期需要分だけ在庫が減少する． 初期在庫（辞書I0）を与えたとき，各期の在庫量が上限（辞書UB）以下，下限（辞書LB)以上でなければいけないとしたとき，段取り費用の合計を最小にする生産計画をたてよ．

S={"0":0,"a":30,"b":50} #S[P,T]：単位生産量　
UB={"0":0,"a":50,"b":50} #UB[p,t]：在庫量の上限 
LB={"0":0,"a":10,"b":10}  #LB[p]：在庫量の下限
I0={"0":0,"a":10,"b":30} #I0[p]:初期在庫

#D[p,t]：需要量
D={('0',1):0,('0',2):0,('0',3):0,('0',4):0,('0',5):0,
   ('a',1):10,('a',2):10,('a',3):30,('a',4):10,('a',5):10,
   ('b',1):20,('b',2):10,('b',3):20,('b',4):10,('b',5):10}

#F[p,q]: 製品p,q間の段取り費用
F={('0',"a"):10,('0',"b"):10,
   ('a',"0"):10,('a',"b"):30,
   ('b',"0"):10,('b',"a"):10}
#
prod=["0","a","b"]      #製品の種類
T=5                   #計画期間は5期

S={"0":0,"a":30,"b":50} #S[P,T]：単位生産量　
UB={"0":0,"a":50,"b":50} #UB[p,t]：在庫量の上限 
LB={"0":0,"a":10,"b":10}  #LB[p]：在庫量の下限
I0={"0":0,"a":10,"b":30} #I0[p]:初期在庫

#D[p,t]：需要量
D={('0',1):0,('0',2):0,('0',3):0,('0',4):0,('0',5):0,
   ('a',1):10,('a',2):10,('a',3):30,('a',4):10,('a',5):10,
   ('b',1):20,('b',2):10,('b',3):20,('b',4):10,('b',5):10}

#F[p,q]: 製品p,q間の段取り費用
F={('0',"a"):10,('0',"b"):10,
   ('a',"0"):10,('a',"b"):30,
   ('b',"0"):10,('b',"a"):10}

model=Model()

X={}          #X[p,t]：製品pを期tに生産するかどうかの0-1変数           
for t in range(1,T+1):
    X[t]=model.addVariable("X{0}".format(t),prod)
    
#constraint                    
for p in prod:
    if p=="0":
        pass
    else:
        for t in range(1,T+1):
            D_temp=0
            for i in range(1,t+1):
                D_temp+=D[p,i]
            con1=Linear("LB{0}_{1}".format(p,t),"inf",LB[p]-I0[p]+D_temp,">=")
            for i in range(1,t+1):
                con1.addTerms(S[p],X[i],p)
            model.addConstraint(con1)

for p in prod:
    if p=="0":
        pass
    else:
        for t in range(1,T+1):
            D_temp=0
            for i in range(1,t+1):
                D_temp+=D[p,i]
            con2=Linear("UB{0}_{1}".format(p,t),"inf",UB[p]-I0[p]+D_temp,"<=")
            for i in range(1,t+1):
                con2.addTerms(S[p],X[i],p)
            model.addConstraint(con2)
        
for p in prod:
    if p=="0":
        pass
    else:
        for q in prod:
            if q=="0" or p==q:
                pass
            else:
                for t in range(2,T+1):
                    con3=Quadratic("obj{0}_{1}_{2}".format(p,q,t),1,0,"<=")
                    con3.addTerms(F[p,q],X[t-1],p,X[t],q)
                    model.addConstraint(con3)


model.Params.TimeLimit=1
sol,violated=model.optimize()

if model.Status==0:
    print ("solution")
    for x in sol:
        print (x,sol[x])
    print ("violated constraint(s)")
    for v in violated:
        print (v,violated[v])

 ================ Now solving the problem ================ 

solution
X1 a
X2 b
X3 a
X4 a
X5 0
violated constraint(s)
obja_b_2 30
objb_a_3 10