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

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

SCOPのトライアルバージョンは， http://logopt.com/scop2/ からダウンロードもしくはgithub https://github.com/scmopt/scop からクローンできる．また，テクニカルドキュメントは， https://scmopt.github.io/manual/14scop.html にある．

重み付き制約充足問題

ここでは，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では数（制約）名や値を文字列で区別するため重複した名前を付けることはできない．なお，使用できる文字列は, 英文字 (a--z, A--Z), 数字 (0--9), 大括弧 ([ ]), アンダーバー (_), および @ に限定される．

それ以外の文字列はすべてアンダーバー (_)に置き換えられる．

Parametersクラス

Parametersクラスで設定可能なパラメータは，以下の通り．

TimeLimit は制限時間を表す．制限時間は正数値を設定する必要があり，その既定値は 600秒である．
OutputFlag は出力フラグを表し，最適化の過程を出力する際の詳細さを制御するためのパラメータである．真(True もしくは正の値)に設定すると詳細な情報を出力し，偽(False もしくは 0)に設定すると最小限の情報を出力する．既定値は偽(0)である．
RandomSeed は乱数の種である．SCOPでは探索にランダム性を加味しているので，乱数の種を変えると，得られる解が変わる可能性がある．乱数の種の既定値は 1である．
Target は制約の逸脱量が目標値以下になったら自動終了させるためのパラメータである．既定値は 0 である．
Initial は，前回最適化の探索を行った際の最良解を初期値とした探索を行うとき True ，それ以外のとき False を表すパラメータである．既定値は False である．最良解の情報は，「変数:値」を1行としたテキストとしてファイル名 scop_best_data.txt に保管されている．このファイルを書き換えることによって，異なる初期解から探索を行うことも可能である．

  optimization will terminate if the solver determines that the optimum penalty value
  for the model is worse than the specified "Target." Non-negative integer. Default = 0.

Parametersクラスの使用例

params = Parameters()
params.TimeLimit = 3 
print(params)

 TimeLimit = 3 
 OutputFlag = 0 
 RandomSeed = 1 
 Taeget = 0 
 Initial = False

変数クラス Variable

変数クラス Variable のインスタンスは，モデルインスタンスの addVariable もしくは addVariables メソッドを用いて生成される．

変数インスタンス=model.addVariable(name, domain）

引数の name は変数名を表す文字列であり， domain は領域を表すリストである．

変数インスタンスのリスト=model.addVariables(names, domain）

引数の names は変数名を要素としたリストであり， domain は領域を表すリストである．

変数クラスは，以下の属性をもつ．

name は変数の名称である．
domain は変数の領域(domain)を表すリストである．変数には領域に含まれる値(value)のうちの1つが割り当てられる．
value は最適化によって変数に割り当てられた値である．最適化が行われる前には None が代入されている．

また，変数インスタンスは，変数の情報を文字列として返すことができる．

Variableクラスの使用例

var = Variable("X[1]", [1,2,3])
print(var)

variable X[1]:['1', '2', '3'] = None

var1 = Variable(domain = [1,2,3])
var2 = Variable(domain = [4,5,6])
print(var1)
print(var2)

variable __x0:['1', '2', '3'] = None
variable __x1:['4', '5', '6'] = None

try:
    var1 = Variable(name = 1, domain = [1,2,3])
except ValueError as error:
    print(error)

Variable name must be a string

モデルクラス

PythonからSCOPをよび出して使うときに，最初にすべきことはモデルクラス Model のインスタンスを生成することである．たとえば， 'test' と名付けたモデルインスタンス model を生成したいときには，以下のように記述する．

from scop import *
model = Model('test')

インスタンスの引数はモデル名であり，省略すると無名のモデルが生成される．

Model クラスは，以下のメソッドをもつ．

addVariable(name，domain) はモデルに1つの変数を追加する．引数の name は変数名を表す文字列であり，domain は領域を表すリストである．変数名を省略すると，自動的に __x[ 通し番号 ] という名前が付けられる．

領域を定義するためのリスト domain の要素は，文字列でも数値でもかまわない．（ただし内部表現は文字列であるので，両者は区別されない．）

以下に例を示す．

x = model.addVarriable('var')                     # domain  is set to [] 
x = model.addVariable(name='var',domain=[1,2,3])  # arguments by name
x = model.addVariable('var',['A','B','C'])        # arguments by position

1行目の例では，変数名を 'var' と設定した空の領域をもつ変数を追加している． 2行目の例では，名前付き引数で変数名と領域を設定している．領域は 1,2,3 の数値である． 3行目の例では，領域を文字列として変数を追加している．

addVariables(names,domain) はモデルに，同一の領域をもつ複数の変数を同時に追加する．引数の names は変数名を要素としたリストであり，domain は領域を表すリストである．領域を定義するためのリストの要素は，文字列でも数値でもかまわない．
addConstriant(con) は制約インスタンス con をモデルに追加する．制約インスタンスは，制約クラスを用いて生成されたインスタンスである．制約インスタンスの生成法については，以下で解説する．
optimize はモデルの求解（最適化）を行うメソッドである．最適化のためのパラメータは，パラメータ属性 Params で設定する．返値は，最適解の情報を保管した辞書と，破った制約の情報を保管した辞書のタプルである．

たとえば，以下のプログラムでは最適解を辞書 sol に，破った制約を辞書 violated に保管する．

sol,violated= model.optimize()

最適解や破った制約は，変数や制約の名前をキーとし，解の値や制約逸脱量を値とした辞書であるので，最適解の値と逸脱量を出力するには，以下のように記述すれば良い．

for x in sol:
    print(x, sol[x])
for v in violated:
    print(v,violated[v])

　引数：

　- cloud: 複数人が同時実行する可能性があるときTrue（既定値はFalse）; Trueのとき，ソルバー呼び出し時に生成されるファイルにタイムスタンプを追加し，計算終了後にファイルを消去する．

モデルインスタンスは，以下の属性をもつ．

name はモデルの名前である．コンストラクタの引数として与えられる．省略可で既定値は ' ' である．
variables は変数インスタンスのリストである．
constraints は制約インスタンスのリストである．
varDict は制約名をキーとし，変数インスタンスを値とした辞書である．
Params は求解（最適化）の際に用いるパラメータを表す属性を保管する．
Status は最適化の状態を表す整数である．状態の種類と意味を以下の表に示す．

最適化の状態を表す整数と意味

状態の定数	説明
0	最適化成功
1	求解中にユーザが Ctrl-C を入力したことによって強制終了した．
2	入力データファイルの読み込みに失敗した．
3	初期解ファイルの読み込みに失敗した．
4	ログファイルの書き込みに失敗した．
5	入力データの書式にエラーがある．
6	メモリの確保に失敗した．
7	実行ファイル scop.exe のよび出しに失敗した．
10	モデルの入力は完了しているが，まだ最適化されていない．
負の値	その他のエラー

また，モデルインスタンスは，モデルの情報を文字列として返すことができる．

Modelクラスの使用例

model = Model('test')
model.addVariable(name="x[1]",domain=[0,1,2])
print(model)
model.Params.TimeLimit= 1
sol, violated = model.optimize()
print("solution=", sol)
print("violated constraints=", violated)
model.Status

Model:test 
number of variables = 1  
number of constraints= 0  
variable x[1]:['0', '1', '2'] = None

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

solution= {'x[1]': '0'}
violated constraints= {}

0

制約クラス

制約クラスは，以下で定義する線形制約クラス，2次制約クラス，相異制約クラスの基底クラスである．

線形制約クラス

最も基本的な制約は，線形制約である．線形制約は $$ 線形項1 + 線形項2 + \cdots 制約の方向 (\leq, \geq, =) 右辺定数 $$ の形で与えられる．線形項は，「変数」が「値」をとったとき $1$，それ以外のとき $0$ を表す値変数(value variable) $x[変数,値]$ を用いて， $$ 係数 \times x[変数,値] $$ で与えられる．

ここで，係数ならびに右辺定数は整数とし，制約の方向は以下($\leq$)，以上($\geq$)，等しい($=$)から1つを選ぶ必要がある．

線形制約クラス Linear のインスタンスは，以下のように生成する．

線形制約インスタンス=Linear(name, weight=1, rhs=0, direction='<=')

引数の意味は以下の通り．

name は制約の名前を表す．これは制約を区別するための名称であり，固有の名前を文字列で入力する必要がある．（名前が重複した場合には，前に定義した制約が無視される．）名前を省略した場合には，自動的に __CON[ 通し番号 ] という名前が付けられる．これは，以下の2次制約や相異制約でも同じである．
weight は制約の重みを表す．重みは，制約の重要性を表す正数もしくは文字列 'inf' である．ここで 'inf' は無限大を表し，絶対制約を定義するときに用いられる．重みは省略することができ，その場合の既定値は 1である．
rhs は制約の右辺定数(right hand side)を表す．右辺定数は，制約の右辺を表す定数（整数値）である．右辺定数は省略することができ，その場合の既定値は 0 である．
direction は制約の向きを表す．制約の向きは， '<=' , '>=' , '=' のいずれかの文字列とする．既定値は '<=' である．

上の引数はすべて Linear クラスの属性となる．最適化を行った後では，制約の左辺の評価値が属性として参照可能になる．

lhs は制約の左辺の値を表す．これは最適化によって得られた変数の値を左辺に代入したときの評価値である．最適化を行う前には 0 が代入されている

Linear クラスは，以下のメソッドをもつ．

addTerms(coeffs, vars，values) は，線形制約の左辺に1つもしくは複数の項を追加する．

addTerms メソッドの引数の意味は以下の通り．
- coeffs は追加する項の係数もしくは係数リスト．係数もしくはリストの要素は整数．
- vars は追加する項の変数インスタンスもしくはの変数インスタンスのリスト．リストの場合には，リストcoeffsと同じ長さをもつ必要がある．
- values は追加する項の値もしは値のリスト．リストの場合には，リストcoeffsと同じ長さをもつ必要がある．　
addTerms メソッドは，1つの項を追加するか，複数の項を一度に追加する．1つの項を追加する場合には，引数の係数は整数値，変数は変数インスタンスで与え，値は変数の領域の要素とする．複数の項を一度に追加する場合には，同じ長さをもつ，係数，変数インスタンス，値のリストで与える．

たとえば，項をもたない線形制約インスタンス L に対して，

L.addTerms(1, y, 'A')

と1つの項を追加すると，制約の左辺は

1 x[y, 'A']

となる．ここで x は値変数（ y が 'A' になるとき 1，それ以外のとき 0 の仮想の変数）を表す．

同様に，項をもたない線形制約インスタンス L に対して，

L.addTerms([2, 3, 1], [y, y, z], ['C', 'D', 'C'])

と3つの項を同時に追加すると，制約の左辺は以下のようになる．

 2 x[y,'C']  + 3 x[y,'D'] + 1 x[z,'C']

setRhs(rhs) は線形制約の右辺定数を rhs に設定する．引数は整数値であり，既定値は 0 である．
setDirection(dir) は制約の向きを設定する．引数 dir は '<=' , '>=' , '=' のいずれかの文字列とする．既定値は '<=' である．
setWeight(weight) は制約の重みを weight に設定する．引数は正数値もしくは文字列 'inf' である．ここで 'inf' は無限大を表し，絶対制約を定義するときに用いられる．

また，線形制約クラス Linear は，制約の情報を文字列として返すことができる．

Linearクラスの使用例

L = Linear(name = "a linear constraint", rhs = 10, direction = "<=" )
x = Variable(name="x", domain=["A","B","C"] )
L.addTerms(3, x, "A")
print(L)

a_linear_constraint: weight= 1 type=linear 3(x,A) <=10

try:
    L = Linear(name = "a linear constraint", rhs = 10.56, direction = "=" )
except ValueError as error:
    print(error)
    
x = Variable(name="x", domain=["A","B","C"] )
try:
    L.addTerms(3.1415, x, "A")
except ValueError as error:
    print(error)

Right-hand-side must be an integer.
Coefficient must be an integer.

2次制約クラス

SCOPでは（非凸の）2次関数を左辺にもつ制約（2次制約）も扱うことができる．

2次制約は， $$ 2次項1 + 2次項2 + \cdots 制約の方向 (\leq, \geq, =) 右辺定数 $$ の形で与えられる．ここで2次項は， $$ 係数 \times x[変数1,値1] \times x[変数2,値2] $$ で与えられる．

2次制約クラス Quadratic のインスタンスは，以下のように生成する．

2次制約インスタンス=Quadratic(name, weight=1, rhs=0, direction='<=')

2次制約クラスの引数と属性は，線形制約クラス Linear と同じである．

Quadratic クラスは，以下のメソッドをもつ．

addTerms(coeffs,vars1,values1,vars2,values2) は2次制約の左辺に2つの変数の積から成る項を追加する．

2次制約に対する addTerms メソッドの引数は以下の通り．

coeffs は追加する項の係数もしくは係数のリスト．係数もしくはリストの要素は整数．
vars1 は追加する項の第1変数インスタンスもしくは変数インスタンスのリスト．リストの場合には，リスト coeffs と同じ長さをもつ必要がある．
values1 は追加する項の第1変数の値もしくは値のリスト．リストの場合には，リスト coeffs と同じ長さをもつ必要がある．
vars2 は追加する項の第2変数の変数インスタンスもしくは変数インスタンスのリスト．リストの場合には，リスト coeffs と同じ長さをもつ必要がある．
values2 は追加する項の第2変数の値もしくは値のリスト．リストの場合には，リスト coeffs と同じ長さをもつ必要がある．

addTerms メソッドは，1つの項を追加するか，複数の項を一度に追加する． 1つの項を追加する場合には，引数の係数は整数値，変数は変数インスタンスで与え，値は変数の領域の要素とする．複数の項を一度に追加する場合には，同じ長さをもつ，係数，変数インスタンス，値のリストで与える．

たとえば，項をもたない2次制約インスタンス Q に対して，

Q.addTerms(1, y, 'A', z, 'B')

と1つの項を追加すると，制約の左辺は

1 x[y,'A'] * x[z,'B']

となる．

同様に，項をもたない2次制約インスタンス Q に対して，

Q.addTerms([2, 3, 1], [y, y, z], ['C', 'D', 'C'], [x, x, y], ['A', 'B', 'C'])

と3つの項を同時に追加すると，制約の左辺は以下のようになる．

2 x[y,'C'] * x[x,'A'] + 3 x[y,'D'] * x[x,'B'] + 1 x[z,'C'] * x[y,'C']

setRhs(rhs) は2次制約の右辺定数を rhs に設定する．引数は整数値であり，既定値は 0 である．
setDirection(dir) は制約の向きを設定する．引数 dir は '<=' , '>=' , '=' のいずれかの文字列とする．既定値は '<=' である．
setWeight(weight) は制約の重みを設定する．引数は正数値もしくは文字列 'inf' である． 'inf' は無限大を表し，絶対制約を定義するときに用いられる．

また，2次制約クラス Quadratic は，制約の情報を文字列として返すことができる．

Quadraticクラスの使用例

Q = Quadratic(name = "a quadratic constraint", rhs = 10, direction = "<=" )
x = Variable(name="x",domain=["A","B","C"] )
y = Variable(name="y",domain=["A","B","C"] )
Q.addTerms([3,9], [x,x], ["A","B"], [y,y], ["B","C"])
print(Q)

a_quadratic_constraint: weight=1 type=quadratic 3(x,A)(y,B) 9(x,B)(y,C) <=10

相異制約クラス

相異制約は，変数の集合に対し, 集合に含まれる変数すべてが異なる値をとらなくてはならないことを規定する．これは組合せ的な構造に対する制約であり，制約最適化の特徴的な制約である．

SCOPにおいては，値が同一であるかどうかは，値の名称ではなく，変数のとりえる値の集合（領域）を表したリストにおける順番（インデックス）によって決定される. たとえば, 変数 var1 および var2 の領域がそれぞれ ['A','B'] ならびに ['B','A'] であったとき，変数 var1 の値 'A', 'B' の順番はそれぞれ 0 と 1，変数 var2 の値 'A', 'B' の順番はそれぞれ 1 と 0 となる．したがって，相異制約を用いる際には変数に同じ領域を与えることが（混乱を避けるという意味で）推奨される．

相異制約クラス Alldiff のインスタンスは，以下のように生成する．

相異制約インスタンス = Alldiff(name, varlist, weight)

引数の名前と既定値は以下の通り．

name は制約名を与える．
varlist は相異制約に含まれる変数インスタンスのリストを与える．これは，値の順番が異なることを要求される変数のリストであり，省略も可能である．その場合の既定値は，空のリストとなる．ここで追加する変数は，モデルクラスに追加された変数である必要がある．
weight は制約の重みを与える．

相異制約の制約名と重みについては，線形制約クラス Linear と同じように設定する．上の引数は Alldiff クラスの属性でもある．その他の属性として最適化した後で得られる式の評価値がある．

lhs は左辺(left hand side)の評価値を表し，最適化された後に，同じ値の番号（インデックス）をもつ変数の数が代入される．

Alldiff クラスは，以下のメソッドをもつ．

addVariable(var) は相異制約に1つの変数インスタンス var を追加する．
addVariables(varlist) は相異制約の変数インスタンスを複数同時に（リスト varlist として）追加する．
setWeight(weight) は制約の重みを設定する．引数は正数値もしくは文字列 'inf' である． 'inf' は無限大を表し，絶対制約を定義するときに用いられる．

また，相異制約クラス Alldiff は，制約の情報を文字列として返すことができる．

Alldiffクラスの使用例

A = Alldiff(name="a alldiff constraint")
x = Variable(name="x",domain=["A","B","C"] )
y = Variable(name="y",domain=["A","B","C"] )
A.addVariables([x,y])
print(A)

a_alldiff_constraint: weight= 1 type=alldiff  y x ;

最適化の描画関数 plot_scop

SCOPはメタヒューリスティクスによって解の探索を行う．一般には，解の良さと計算時間はトレードオフ関係がある．つまり，計算時間をかければかけるほど良い解を得られる可能性が高まる．どの程度の計算時間をかければ良いかは，最適化したい問題例（問題に数値を入れたもの）による． plot_scopは，横軸に計算時間，縦軸に目的関数値をプロットする関数であり，最適化を行ったあとに呼び出す．得られるPlotlyの図は，どのくらいの計算時間をかければ良いかをユーザーが決めるための目安を与える．

たとえば以下の例の図から，500秒程度の計算時間で良好な解を得ることができるが，念入りにさらに良い解を探索したい場合には2000秒以上の計算時間が必要なことが分かる．

fig = plot_scop()
#plotly.offline.plot(fig);

例題

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

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

from scop import *

例題仕事の割当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  B C A ;  :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  C A 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$ の問題を以下の手順で問題を解く．

各マス目 $(i,j), i=0,1,2, j=0,1,2$ に対して変数 $x[i,j]$ を準備して，その領域を $1$ から $9$ までの数とする．
各マス目には異なる数字を入れる必要があるので，すべての変数のリストを入れた相異制約 (Alldiff) を追加する．この制約は絶対制約とする．
さらに，各行（$i=0,1,2$)と各列($j=0,1,2$)に対して，その和がちょうど $15 = (1+2+\cdots+9)/3$ になるという制約を追加する．これらの制約は考慮制約とし，逸脱の重みは $1$ とする．
最適化を行い，解を表示する．

行の集合を $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]]

デバッグのコツ

ここでは，デバッグのコツについて述べる．

SCOPでは，変数や制約のオブジェクトをprint関数で表示することができる．プログラムを少しずつ書き， printでオブジェクトが意図通りになっているかを出力して確認する．

SCOPは，Pythonからテキストファイルを生成して，それを実行ファイルに入れて求解する．入力ファイル名は scop_input.txt で出力ファイル名は scop_output.txt である．これらのファイルをエディタでみて，意図通りかどうかを確認することもできる．

その際に，大きな問題例だとデバッグしにくいので，非常に小さい例題を準備することを推奨する．また，制約条件も簡単なものから順に複雑な制約に変更していき，実行不能になったときの原因が容易に分かるようにすべきである．

大規模な実際問題の解決法

大きな問題例(問題に数値を与えたもの；instance)を直接モデル化して実験するのではなく，小規模ですべての機能をテストできるような問題例を準備すべき．
SCOPでは，データはすべて整数値にしなければならない．
1つの式の評価が大きな整数になることを避ける．
大きな整数になる可能性がある項の和に対する線形制約のかわりに，項を分解してそれぞれの式を評価する．
以前探索を行ったときの最良解の情報は，「変数:値」を1行としたテキストとしてファイル名 scop_best_data.txt に保管されている．これを初期解として再最適化を行うには， model.Params.InitialをTrueに設定する．これによって2回め以降の探索の高速化が可能になる．また，このファイルを書き換えることによって，異なる初期解から探索を行うこともできる．

ターミナル（コマンドプロンプト）での実行

SCOPは，Pythonからテキストファイルを生成して，それを実行ファイルに入れて求解する．入力ファイル名は scop_input.txt で出力ファイル名は scop_output.txt である．デバッグ時には，入力ならびに出力ファイルを適当なテキストエディタで開いてみると良い．

また，デバッグ時だけでなく，大規模問題例を解くときや，アルゴリズムのログを確認しながら実験したい場合には，ターミナル（コマンドプロンプト）からの実行を推奨する．

コマンドはOSによって異なるが，scop(Mac), scop-linux (linux), scop.exe (Windows) を以下の形式で実行する．

./scop       < scop_input.txt

./scop-linux < scop_input.txt

scop         < scop_input.txt

オプションの引数は，「-オプション名引数」の形式で追加する． scop --help で表示されるオプションは，以下の通りである．

  -noadjust         deactivate weight adjustment mechanism
  -display #        set log display level
  -iteration #      set iteration limit
  -initsolfile file set initial solution file name
  -inputfile file   set input file name
  -interval #       set log display interval
  -logfile file     set log file name
  -noimprovement #  set iteration limit for no improvement
  -seed #           set random seed
  -logfile file     set log file name
  -target #         set target
  -time #.#         set cpu time limit in second

たとえば，実行時のログを1反復ごとに表示させたい場合には，

./scop     -interval 1 < scop_input.txt

./scop-linux -interval 1 < scop_input.txt

scop         -interval 1 < scop_input.txt

とすれば良い．

練習問題

問題多制約ナップサック

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

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

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

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

問題最大安定集合

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

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

adj=[[2],[3],[0,3,4,5],[1,2,5],[2],[2,3]]

問題グラフ彩色

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

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

adj=[[2],[3],[0,3,4,5],[1,2,5],[2],[2,3]]

問題グラフ分割

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

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

adj=[[1,4],[0,2,4],[1],[4,5],[0,1,3,5],[3,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]]

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

問題ビンパッキング

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

問題最適化版の$8$-クイーン

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

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

問題 $k$-メディアン

顧客から最も近い施設への距離の「合計」を最小にするようにグラフ内の点上または空間内の任意の点から施設を選択する問題である．メディアン問題においては，選択される施設の数があらかじめ決められていることが多く，その場合には選択する施設数 $k$ を頭につけて $k$-メディアン問題($k$-median problem)とよばれる．

顧客数を $n$ とし，顧客の集合を $I$，施設の配置可能な点の集合を $J$ とする．簡単のため $I=J$ とし， 2地点間の移動距離を以下のデータ生成関数によって生成したものとしたとき，地点数 $|I|=|J|=200$， $k=20$ の $k$-メディアン問題の解を求めよ．

注意：この問題はトライアル版では解けない．

import random
import math
import networkx as nx

def distance(x1,y1,x2,y2):
    return math.sqrt((x2-x1)**2 + (y2-y1)**2)

random.seed(67)
def make_data(n,m,same=True):
    if same == True:
        I = range(n)
        J = range(m)
        x = [random.random() for i in range(max(m,n))]  
        # positions of the points in the plane
        y = [random.random() for i in range(max(m,n))]
    else:
        I = range(n)
        J = range(n,n+m)
        x = [random.random() for i in range(n+m)]    
        # positions of the points in the plane
        y = [random.random() for i in range(n+m)]
    c = {}
    for i in I:
        for j in J:
            c[i,j] = distance(x[i],y[i],x[j],y[j])

    return I,J,c,x,y

n = 200
m = n
I, J, c, x_pos, y_pos = make_data(n,m,same=True)
k = 20

問題数独パズル

与えられる盤面は 3 × 3 のブロックに区切られた 9 × 9 の格子である.その幾つかのマスには，あらかじめ 1 から 9 までの数字が入っている. 与えられた盤面の空いているマスに次の規則を満たすように 1 から 9 までの数字を入れる.

どの横行，縦列，および太線で囲まれた 3 × 3 のブロックにも同じ数字が入ってはいけない.

これは，そのどれらにもちょうど 1 つずつ 1 から 9 までの数字が入ることと同じ意味である.

注意：この問題はトライアル版では解けない．

grid1 = """
.14.5...3
6....942.
8..1...9.
..5.9..4.
4..7.8..2
.7..2.6..
.9...1..5
.283....4
5...6.71.
"""
def values_from_grid(grid):
    values = []
    digits = "123456789"
    chars = [c for c in grid if c in digits or c in '0.']
    grid_int = map(lambda x: int(x) if x != "." else 0, chars)

    count = 0
    row = []
    for i in grid_int:
        row.append(i)
        count += 1
        if count % 9 == 0: #行毎に分割
            values.append(row)
            row = []
    return values
value = values_from_grid(grid1)
value

[[0, 1, 4, 0, 5, 0, 0, 0, 3],
 [6, 0, 0, 0, 0, 9, 4, 2, 0],
 [8, 0, 0, 1, 0, 0, 0, 9, 0],
 [0, 0, 5, 0, 9, 0, 0, 4, 0],
 [4, 0, 0, 7, 0, 8, 0, 0, 2],
 [0, 7, 0, 0, 2, 0, 6, 0, 0],
 [0, 9, 0, 0, 0, 1, 0, 0, 5],
 [0, 2, 8, 3, 0, 0, 0, 0, 4],
 [5, 0, 0, 0, 6, 0, 7, 1, 0]]

問題シフトスケジューリング

あなたは，24時間営業のハンバーガーショップのオーナーであり，スタッフの1週間のシフトを組むことに頭を悩ませている．スタッフの時給は同じであると仮定したとき，以下の制約を満たすシフトを求めよ．

1シフトは 8時間で，朝，昼，晩の3シフトの交代制とする．
3人のスタッフは，1日に高々1つのシフトしか行うことができない．
繰り返し行われる1週間のスケジュールの中で，スタッフは最低3日間は勤務しなければならない．
スタッフの夜勤は最大で4回とする．
各スタッフは1日以上休みを取る必要がある．
各シフトに割り当てられるスタッフの数は，ちょうど 1人でなければならない．
異なるシフトを翌日に行ってはいけない．（すなわち異なるシフトに移るときには，必ず休日を入れる必要がある．）
シフト2, 3は，少なくとも2日間は連続で行わなければならない．

注意：この問題はトライアル版では解けない．

問題車の投入順決定

コンベア上に一直線に並んだ車の生産ラインを考える．このラインは，幾つかの作業場から構成され，それぞれの作業場では異なる作業が行われる．いま，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$つ入っているという制約を課すことに相当する

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

１つの生産ラインで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}

制約最適化システム SCOP

重み付き制約充足問題

SCOPの基本クラス

注意

Parametersクラス

class Parameters[source]

Parametersクラスの使用例

変数クラス Variable

class Variable[source]

Variableクラスの使用例

モデルクラス

class Model[source]

Modelクラスの使用例

制約クラス

class Constraint[source]

線形制約クラス

class Linear[source]

Linearクラスの使用例

2次制約クラス

class Quadratic[source]

Quadraticクラスの使用例

相異制約クラス

class Alldiff[source]

Alldiffクラスの使用例

最適化の描画関数 plot_scop

plot_scop[source]

例題

例題 仕事の割当1

例題 仕事の割当2

例題 仕事の割当3

例題 魔方陣

デバッグのコツ

大規模な実際問題の解決法

ターミナル（コマンドプロンプト）での実行

練習問題

問題 多制約ナップサック

問題 最大安定集合

問題 グラフ彩色

問題 グラフ分割

問題 巡回セールスマン

問題 ビンパッキング

問題 最適化版の$8$-クイーン

問題 2次割当

問題 $k$-メディアン

問題 数独パズル

問題 シフトスケジューリング

問題 車の投入順決定

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

`class` `Parameters`[source]

`class` `Variable`[source]

`class` `Model`[source]

`class` `Constraint`[source]

`class` `Linear`[source]

`class` `Quadratic`[source]

`class` `Alldiff`[source]

`plot_scop`[source]

例題仕事の割当1

例題仕事の割当2

例題仕事の割当3

例題魔方陣

問題多制約ナップサック

問題最大安定集合

問題グラフ彩色

問題グラフ分割

問題巡回セールスマン

問題ビンパッキング

問題最適化版の$8$-クイーン

問題数独パズル

問題シフトスケジューリング

問題車の投入順決定

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