ここで提案する SCML (Suply Chain Modeling Language) とは， 様々なサプライ・チェイン最適化モデルを 統一的に表現するための基本言語である．
サプライ・チェインを対象とした最適化モデルには，様々なものがある． 同じ現実問題を対象としていても，その抽象化のレベルによって，様々なモデルが作成できる． エンドユーザーが使用することを想定した場合には， 実ロジスティクス・オブジェクト（トラックや機械などを指す造語）を用いたモデルが有効であり， 一方，システムを作成するロジスティクス・エンジニアや研究者が使用することを想定した場合には， 抽象ロジスティクス・オブジェクト（資源やネットワークなどを指す造語）を用いたモデルが有効である．
また，同じ抽象モデルでも，グラフ・ネットワークなど高度に抽象化されたものから， 時間枠付き配送計画や資源制約付きスケジューリングなど具体的な応用を対象としたものまで 様々な階層に分けることができ，それぞれ利点・弱点がある．
抽象度の階層の最下層にあるのが，実際の現場で用いられるシステムに内在するモデルであり， それらは現実モデルに直結している．最下層の現実モデルの上位にあるのが，様々な典型的なサプライ・チェイン最適化モデル であり，個々のモデルに対して多くの研究がなされている． この階層のモデルの例として，ロジスティクス・ネットワーク設計モデル，安全在庫配置モデル，スケジューリングモデル，配送計画モデルなどがあげられる． ここで考えるSCMLは，これらのモデルを統一的に表現するために設計されている．
サプライ・チェインの理論や実務については，様々な場所で分野を超えた交流が成されているが， 分野ごとに使用されている用語の意味や定義が異なるため，議論が噛み合わず，有益な交流ができているとは言い難い． たとえば，ある分野ではサプライ・チェイン最適化とは「在庫」を適正化するための手法として捉えられており， 別の分野ではサプライ・チェイン最適化とは工場内のスケジューリングを指していたりする．
異分野間の交流は，サプライ・チェインのような複合的な学問体系にとっては必要不可欠なものであるが， これらの現象を目のあたりにし，研究者と実務家が同じ土俵で議論するためのモデルの必要性を痛切に感じた． これが，SCMLを考えるに至った動機である．
ここで提案する抽象モデルでは，サプライ・チェインを， 空間（点と枝からなるネットワーク）， 時間（期），製品，資源，活動（とそれを実行する方法であるモード）などの基本構成要素 から構成されるものと捉える． 構成要素の中心に活動があり， 活動（とそれに付随するモード）が資源を用いて製品を時・空間内に移動させることが サプライ・チェインの本質であり，目的であると考える．
以下では，これらの基本構成要素 (entity) をPython/Pydanticを用いたクラスとして表現していく． また，これらを組み込んだモデルクラスを設計し，最適化や可視化のためのメソッドを準備する．
準備のために基本となるEntityクラスを作っておく。PydanticのBaseModelから派生させる。
Entity間に，集約（グループ化）と非集約関係を定義するために，集約したEntityを表す親 (parent) と，非集約したときのEntityのリストを表す子リスト (children) を定義しておく。 たとえば，月が親のとき子はその月に含まれる週や日になる。
属性
親子関係の定義には，makeGroup関数などを用いる。
Entity (name:Union[int,str], parent:Union[int,str,NoneType]=None, children:Optional[List[Union[int,str]]]=None)
entity1 = Entity(name="entity1")
entity2 = Entity(name="entity2")
entity3 = Entity(name="entity3")
entity1.setParent(entity2)
print(entity2)
entity1.addChildren([entity2, entity3])
print(entity1)name='entity2' parent=None children=['entity1']
name='entity1' parent='entity2' children=['entity2', 'entity3']期 (Period) は，時間をモデル化するために用いられる連続時間を離散化したものである． Entityクラスから派生させるので，name, parent, chidren属性をもつ（以下同様）。
最も単純な期集合は，有限な正数 \(T\)，時刻の列 \(0= t_1 < t_2 < \cdots < t_T\) を与えたとき， 区間 \((t_i,t_{i+1}] ~(i=1,\cdots,T-1)\) の順序付き集合として生成される． \(t_{i+1}-t_i\) がすべて同じであるとき，\(t_{i+1}-t_i (=\delta)\) を期集合の幅とよぶ． サプライ・チェインをモデル化する際には， 意思決定レベルの違いにより様々な幅をもつ期集合が定義される． ここでは，それらを集めたものを「期」と定義する． 期（集合）に対しても，点と同様に集約・非集約関係が定義できる． たとえば，日を集約したものが週であり，週を集約したものが月（もしくは四半期や年）となる．
期は，名前 name と開始時刻を表す start を引数として生成される． startは整数，浮動小数点数，日付時刻型，日付型，時刻型のいずれかの型をもつ． モデルに追加された期は非減少順（同じ場合には名前の順）に並べ替えられ， 小さい順に \(0,1,2,\cdots\) とインデックスが付けられる．
モデルの最適化の際には，開始期と終了期を指定して，一部の期間に対してだけ最適化することができる．
属性
Modelに付加された期インスタンスは，startの非減少順（同点の場合には名前順）に並べられ，\(0\) から始まる整数値のインデックスが付与される． 最適化は，モデルに与えられた最適化期間（省略された場合はすべての期間）に対して行われる．
Period (name:Union[int,str], parent:Union[int,str,NoneType]=None, children:Optional[List[Union[int,str]]]=None, start:Union[int,flo at,datetime.datetime,datetime.date,datetime.time,NoneType])
Resourceは，有限な資源を表すクラスである。 サプライ・チェインを構成する企業体は， 生産ライン，機械，輸送機器（トラック，船，鉄道，飛行機）， 金（財務資源），人（人的資源）などの資源から構成される． 資源集合に対しても，点と同様に集約・非集約関係が定義できる．
属性
メソッド
Resource (name:Union[int,str], parent:Union[int,str,NoneType]=None, children:Optional[List[Union[int,str]]]=None, rtype:Optional[str]=None, capacity:Union[int,float,Dict,NoneType]=None, fixed_cost:Union[int,float,NoneType]=0, modes:Optional[Dict[str,Tuple]]=None)
Productは，製品（品目，半製品，原材料，部品，中間製品，完成品など）を表すクラスである． 製品は，外部から供給され，系内を姿形を変えながら移動した後，外部へ出て行く「もの」 である．製品は，ネットワークフローのモデルの用語では，品種（commodity）とよばれ，資源とは区別される．
属性
Product (name:Union[int,str], parent:Union[int,str,NoneType]=None, children:Optional[List[Union[int,str]]]=None, weight:Optional[float]=0.0, volume:Optional[float]=0.0, value:Optional[float]=0.0)
Product(name='Prod1', parent=None, children=None, weight=0.0, volume=0.0, value=0.0)Mode（モード，方策）は，以下でt定義する Activityクラスに付随し，活動を行うための具体的な方法を表す． 方策 (policy) は動的最適化の枠組みにおいては，システムの状態 (state) を行動 (action) への写像として定義されるが， ここでは活動 (activity) を実行するための方法として捉えるものとする． また，製品がどの原材料や部品から構成されるかを表す部品展開表 (Bill Of Materials: BOM) の情報ももつ。
属性
関連付けられる活動のタイプが生産(make)の場合には，生産量を表す変数 \(x_{m}^p\) が定義され， 段取り(setup)の場合には，生産量 \(x_{m}^p\) と段取りの有無を表す \(0\)-\(1\) 変数 \(y_{m}^p\) が定義される． ここで\(p\) は活動に付随する製品である．
TODO: 一般化としての区分的線形な費用
メソッド
Mode (name:Union[int,str], parent:Union[int,str,NoneType]=None, children:Optional[List[Union[int,str]]]=None, components:Optional[Dict[str,float]]=None, byproducts:Optional[Dict[str,float]]=None, fixed_requirement:Optional[Dict[str,float]]=None, variable_requirement:Optional[Dict[str,float]]=None, fixed_cost:Optional[float]=0.0, variable_cost:Optional[float]=0.0, piecewise_x:Optional[List[float]]=None, piecewise_y:Optional[List[float]]=None, upper_bound:Optional[float]=None, activities:Optional[set]=None, periods:Optional[Dict[Union[int,str],__main__.Period]]=None, duration:Union[int,float,NoneType]=None, service_time:Union[int,float,NoneType]=None, service_level:Optional[float]=0.9)
mode1 = Mode(name = "Mode1")
mode2 = Mode(name = "Mode2")
mode1.addResource(resource=res1, fixed=1)
mode1.addCost(fixed=100., variable=1.)
print(mode1)name='Mode1' parent=None children=None components=None byproducts=None fixed_requirement={'Res1': 1} variable_requirement={'Res1': 0.0} fixed_cost=100.0 variable_cost=1.0 piecewise_x=None piecewise_y=None upper_bound=None activities=None periods=None duration=None service_time=None service_level=0.9サプライ・チェインとは，資源を時・空間内で消費・生成させることであると捉え， 資源を消費・生成する基本となる単位を活動 (activity) とよぶ． Activityはサプライ・チェインに必要な諸活動（作業，タスク，ジョブ）を表すクラスである．
活動集合に対しても，点と同様に集約・非集約関係が定義できる． また，活動は，点もしくは枝上で定義することもできる。 その際には，活動は局所的に定義され，そうでない場合には大域的に定義される。
属性
メソッド
Activity (name:Union[int,str], parent:Union[int,str,NoneType]=None, children:Optional[List[Union[int,str]]]=None, atype:Optional[str]='make', mtype:Optional[str]='select', product:__main__.Product=None, modes:Optional[Dict[str,__main__.Mode]]=None, nodes:Optional[Set]=None, arcs:Optional[Set]=None)
act1 = Activity(name="Act1", product= prod1)
act2 = Activity(name="Act2", product= prod1)
act1.addMode(mode1)
act1.addModes([mode1,mode2])
act2.addMode(mode1)
print(act1)name='Act1' parent=None children=None atype='make' mtype='select' product=Product(name='Prod1', parent=None, children=None, weight=0.0, volume=0.0, value=0.0) modes={'Mode1': Mode(name='Mode1', parent=None, children=None, components=None, byproducts=None, fixed_requirement={'Res1': 1}, variable_requirement={'Res1': 0.0}, fixed_cost=100.0, variable_cost=1.0, piecewise_x=None, piecewise_y=None, upper_bound=None, activities={'Act1', 'Act2'}, periods=None, duration=None, service_time=None, service_level=0.9), 'Mode2': Mode(name='Mode2', parent=None, children=None, components=None, byproducts=None, fixed_requirement=None, variable_requirement=None, fixed_cost=0.0, variable_cost=0.0, piecewise_x=None, piecewise_y=None, upper_bound=None, activities={'Act1'}, periods=None, duration=None, service_time=None, service_level=0.9)} nodes=None arcs=None原料供給地点，工場，倉庫の配置可能地点，顧客（群），作業工程，在庫の一時保管場所など， サプライ・チェインに関わるすべての地点を総称して点 (node) とよぶ． Nodeは点を表すクラスである． 点集合間には集約・非集約関係が定義できる．たとえば，顧客を集約したものが顧客群となる．
属性 - location: 経度・緯度のタプル - location_index: Matricesクラスで定義された行列のインデックス - activities: 点で定義された活動の辞書（名前をキーとし，活動オブジェクトを値とする
メソッド - addActivity: 活動を追加する。 - addActivities: 活動のリストを一度に追加する。
Node (name:Union[int,str], parent:Union[int,str,NoneType]=None, children:Optional[List[Union[int,str]]]=None, location:Union[NoneTy pe,Annotated[Tuple[float,float],Json],Tuple[float,float]]=None, location_index:Optional[int]=None, activities:Optional[Dict[str,__main__.Activity]]=None)
node1 = Node(name="Node1", location=(123.45, 567.89))
node2 = Node(name="Node2")
node1.addActivity(act1)
node1.model_dump_json(exclude_defaults=True)'{"name":"Node1","location":[123.45,567.89],"activities":{"Act1":{"name":"Act1","product":{"name":"Prod1"},"modes":{"Mode1":{"name":"Mode1","fixed_requirement":{"Res1":1.0},"variable_requirement":{"Res1":0.0},"fixed_cost":100.0,"variable_cost":1.0,"activities":["Act1","Act2"]},"Mode2":{"name":"Mode2","activities":["Act1"]}},"nodes":["Node1"]}}}'点の対（2つ組）を枝 (arc) とよぶ． Arcは枝を表すクラスである．枝は，点と点の関係を表し，空間上の移動を表現するために用いられる．
属性
メソッド
Arc (name:Union[int,str], parent:Union[int,str,NoneType]=None, children:Optional[List[Union[int,str]]]=None, source:__main__.Node, sink:__main__.Node, distance:Union[int,float,NoneType]=None, time:Union[int,float,NoneType]=None, activities:Optional[Dict[str,__main__.Activity]]=None)
arc1 = Arc(name="arc1", source=node1, sink=node2)
arc1.addActivity(act1)
arc1.model_dump_json(exclude_defaults=True)'{"name":"arc1","source":{"name":"Node1","location":[123.45,567.89],"activities":{"Act1":{"name":"Act1","product":{"name":"Prod1"},"modes":{"Mode1":{"name":"Mode1","fixed_requirement":{"Res1":1.0},"variable_requirement":{"Res1":0.0},"fixed_cost":100.0,"variable_cost":1.0,"activities":["Act1","Act2"]},"Mode2":{"name":"Mode2","activities":["Act1"]}},"nodes":["Node1"],"arcs":[["Node1","Node2"]]}}},"sink":{"name":"Node2"},"activities":{"Act1":{"name":"Act1","product":{"name":"Prod1"},"modes":{"Mode1":{"name":"Mode1","fixed_requirement":{"Res1":1.0},"variable_requirement":{"Res1":0.0},"fixed_cost":100.0,"variable_cost":1.0,"activities":["Act1","Act2"]},"Mode2":{"name":"Mode2","activities":["Act1"]}},"nodes":["Node1"],"arcs":[["Node1","Node2"]]}}}'Dataは，製品，ノード，期ごとに定義される数値データを保持するクラスである． 名前は必要ないのでEntityクラスから派生させない． モデル上では，点，製品，期をインデックスとして定義することができる． インデックスに依存しないデータに対しては \(*\)（アスタリスク）の文字列を用いて，モデル内に保管される．
属性
Data (dtype:str='demand', amount:Union[int,float,NoneType]=0, std:Optional[float]=0.0, over_penalty:Union[int,float,NoneType]=999999, under_penalty:Union[int,float,NoneType]=999999)
Usage docs: https://docs.pydantic.dev/2.6/concepts/models/
A base class for creating Pydantic models.
Attributes: class_vars: The names of classvars defined on the model. private_attributes: Metadata about the private attributes of the model. signature: The signature for instantiating the model.
__pydantic_complete__: Whether model building is completed, or if there are still undefined fields.
__pydantic_core_schema__: The pydantic-core schema used to build the SchemaValidator and SchemaSerializer.
__pydantic_custom_init__: Whether the model has a custom `__init__` function.
__pydantic_decorators__: Metadata containing the decorators defined on the model.
This replaces `Model.__validators__` and `Model.__root_validators__` from Pydantic V1.
__pydantic_generic_metadata__: Metadata for generic models; contains data used for a similar purpose to
__args__, __origin__, __parameters__ in typing-module generics. May eventually be replaced by these.
__pydantic_parent_namespace__: Parent namespace of the model, used for automatic rebuilding of models.
__pydantic_post_init__: The name of the post-init method for the model, if defined.
__pydantic_root_model__: Whether the model is a `RootModel`.
__pydantic_serializer__: The pydantic-core SchemaSerializer used to dump instances of the model.
__pydantic_validator__: The pydantic-core SchemaValidator used to validate instances of the model.
__pydantic_extra__: An instance attribute with the values of extra fields from validation when
`model_config['extra'] == 'allow'`.
__pydantic_fields_set__: An instance attribute with the names of fields explicitly set.
__pydantic_private__: Instance attribute with the values of private attributes set on the model instance.Constraintは，制約を定義するためのクラスである． スケジューリングモデルに対しては，系全体での資源制約である再生不能資源とも考えられる．
制約は，活動 \(a\) がモード \(m\) を実行するときに \(1\) の \(0\)-\(1\) 変数 \(x_{am}\) に対する以下の線形制約として記述される．
\[ \sum_{a,m} coeff_{am} x_{am} \leq (=, \geq) rhs \]
制約インスタンスは，以下のメソッドをもつ．
setRhs(rhs)は再生不能資源を表す線形制約の右辺定数をrhsに設定する．引数は整数値（負の値も許すことに注意）とする．
setDirection(dir)は再生不能資源を表す制約の種類をdirに設定する． 引数のdirは’<=‘,’>=‘,’=’のいずれかとする．
addTerms(coeffs,vars,values)は，再生不能資源制約の左辺に1つ，もしくは複数の項を追加するメソッドである． 作業がモードで実行されるときに \(1\)， それ以外のとき \(0\) となる変数（値変数）を x[作業,モード]とすると， 追加される項は， $ 係数 x[作業,モード]
$ と記述される． addTermsメソッドの引数は以下の通り．
制約インスタンスは以下の属性をもつ．
Constraint (name:Union[int,str], parent:Union[int,str,NoneType]=None, children:Optional[List[Union[int,str]]]=None, rhs:Union[int,float,NoneType]=0, direction:Optional[str]='<=', terms:Optional[List[Tuple[Union [int,float],__main__.Activity,__main__.Mode]]]=None, weight:Union[int,float,str,NoneType]='inf')
Modelは上の諸Entityを組み合わせたモデルクラスである。データ入力後にupdateメソッドで最適化の準備を行う．
属性
また，上の諸情報をモデルに付加するメソッドをもつ。
最適化メソッドの引数で，最適化モデルの種類を渡す（‘lnd’, ‘optseq’, ‘risk’, ‘ssa’, ‘shift’, … )
Model (name:Union[int,str], parent:Union[int,str,NoneType]=None, children:Optional[List[Union[int,str]]]=None, activities:Optional[Dict[str,__main__.Activity]]=None, modes:Optional[Dict[str,__main__.Mode]]=None, resources:Optional[Dict[str,__main__.Resource]]=None, products:Optional[Dict[str,__main__.Product]]=None, periods:Optional[Dict[str,__main__.Period]]=None, nodes:Optional[Dict[str,__main__.Node]]=None, arcs:Optional[Dict[str,__main__.Arc]]=None, constraints:Optional[Dict[str,__main__.Constraint]]=None, data:Optional[Dict]=None, interest_rate:Optional[float]=None, period_list:Optional[List[Tuple]]=None, period_index:Optional[Dict[str,int]]=None)
| Type | Details | |
|---|---|---|
| data | Any | |
| Returns | None | type: ignore |
sequenceDiagram
participant Client
participant Cloud
loop
Cloud->>Cloud: Optimization
end
Note right of Cloud: Scheduling <br/> Lot-sizing <br/> Logistics Network Design <br/> Shift <br/> Vehicle Routing
Client->>Cloud: Model Instance in JSON
Cloud->>Client: Solution in JSON
model1 = Model(name="model1")
model1.addData(dtype="demand", product=prod1, amount=100.)
model1.addData(dtype="supply", product=prod1, node=node1, period=period1, amount=100.)
#model1.addValue(product=prod1, node=node1, value=12)
model1.addNode(name="Node1")
con1 = model1.addConstraint(name="constraint1")
con1.setRhs(100)
model1.model_dump_json(exclude_none=True)'{"name":"model1","nodes":{"Node1":{"name":"Node1"}},"constraints":{"constraint1":{"name":"constraint1","rhs":100,"direction":"<=","weight":"inf"}},"data":{"demand,Prod1,*,*":{"dtype":"demand","amount":100.0,"std":0.0,"over_penalty":999999,"under_penalty":999999},"supply,Prod1,Node1,Period1":{"dtype":"supply","amount":100.0,"std":0.0,"over_penalty":999999,"under_penalty":999999}}}'model2 = Model(name="model2")
for i,t in enumerate(pd.date_range(start="2024/1/1", end="2024/1/31", freq="w")):
model2.addPeriod(name= f"period({i})", start=t)
model2.update()
print(model2.period_index)
model2.model_dump_json(exclude_none=True){'period(0)': 0, 'period(1)': 1, 'period(2)': 2, 'period(3)': 3}'{"name":"model2","periods":{"period(0)":{"name":"period(0)","start":"2024-01-07T00:00:00"},"period(1)":{"name":"period(1)","start":"2024-01-14T00:00:00"},"period(2)":{"name":"period(2)","start":"2024-01-21T00:00:00"},"period(3)":{"name":"period(3)","start":"2024-01-28T00:00:00"}},"period_list":[["2024-01-07T00:00:00","period(0)"],["2024-01-14T00:00:00","period(1)"],["2024-01-21T00:00:00","period(2)"],["2024-01-28T00:00:00","period(3)"]],"period_index":{"period(0)":0,"period(1)":1,"period(2)":2,"period(3)":3}}'フローの最適化を行う関数
TODO: 結果をクラスの属性に記入する
モデルの種類によって呼ぶ関数を分岐
optimize (model:__main__.Model, start_period_name:str=None, end_period_name:str=None)
モデルの可視化を行う関数
visualize_network (model:__main__.Model)
visualize (model:__main__.Model, rankdir:str='LR', size:float=5)
いま，顧客数を \(n\)，工場数を \(m\) とし， 顧客を \(i=1,2,\ldots,n\)，工場を \(j=1,2,\ldots,m\) と番号で表すものとする． また，顧客の集合を \(I=\{1,2,\ldots,n \}\)，工場の集合を \(J=\{ 1,2,\ldots,m \}\) とする． 顧客 \(i\) の需要量を \(d_i\)，顧客 \(i\) と施設 \(j\) 間に \(1\) 単位の需要が移動するときにかかる 輸送費用を \(c_{ij}\)，工場 \(j\) の容量 \(M_j\) とする． また， \(x_{ij}\) を工場 \(j\) から顧客 \(i\) に輸送される量を表す連続変数する．
上の記号および変数を用いると，輸送問題は以下の線形最適化問題として定式化できる．
\[ \begin{array}{l l l} minimize & \displaystyle\sum_{i \in I} \displaystyle\sum_{j \in J} c_{ij} x_{ij} & \\ s.t. & \displaystyle\sum_{j \in J} x_{ij} =d_i & \forall i \in I \\ & \displaystyle\sum_{i \in I} x_{ij} \leq M_j & \forall j \in J \\ & x_{ij} \geq 0 & \forall i \in I, j \in J \end{array} \]
目的関数は輸送費用の和の最小化であり，最初の制約は需要を満たす条件， 2番目の制約は工場の容量制約である．
import random
from pprint import pprint
random.seed(123)
m = 2
n = 5
model = Model(name="Transportation")
dummy_product = model.addProduct(name="dummy product")
supplier, customer = {}, {}
arc = {}
activity, mode ={}, {}
for i in range(m):
supplier[i] = model.addNode(name=f"Supplier({i})")
model.addData(dtype="supply", product=dummy_product, node=supplier[i], under_penalty= 0, amount= 100 )
for j in range(n):
customer[j] = model.addNode(name=f"Customer({j})")
model.addData(dtype="demand", product=dummy_product, node=customer[j], amount= 20 )
for i in range(m):
for j in range(n):
arc[i,j] = model.addArc(name=f"arc({i},{j})", source=supplier[i], sink=customer[j])
activity[i,j] = model.addActivity(name=f"act({i},{j})", atype = "transport", product = dummy_product )
activity[i,j].addMode ( Mode(name=f"mode({i},{j})", variable_cost = random.randint(5,10) ) )
arc[i,j].addActivity( activity[i,j] )
#pprint(model.model_dump_json(exclude=None))
gp_model = optimize(model)計画期間 = 0 1
Opt. Val= 600.0容量制約付き施設配置問題(capacitated facility location problem)は，以下のように定義される問題である．
顧客数を \(n\)， 施設数を \(m\) とし， 顧客を \(i=1,2,\ldots,n\)， 施設を \(j=1,2,\ldots,m\) と番号で表すものとする． また，顧客の集合を \(I=\{1,2,\ldots,n\}\)，施設の集合を \(J=\{1,2,\ldots,m\}\) と記す．
顧客 \(i\) の需要量を \(d_i\)，顧客 \(i\) と施設 \(j\) 間に \(1\) 単位の需要が移動するときにかかる 輸送費用を \(c_{ij}\)，施設 \(j\) を開設するときにかかる固定費用を \(f_j\)，容量を \(M_j\) とする．
以下に定義される連続変数 \(x_{ij}\) および \(0\)-\(1\) 整数変数 \(y_j\) を用いる．
\[ x_{ij}= 顧客 i の需要が施設 j によって満たされる量 \]
\[ y_j = \left\{ \begin{array}{ll} 1 & 施設 j を開設するとき \\ 0 & それ以外のとき \end{array} \right. \]
上の記号および変数を用いると，容量制約付き施設配置問題は以下の混合整数最適化問題 として定式化できる．
\[ \begin{array}{l l l } minimize & \sum_{j \in J} f_j y_j +\sum_{i \in I} \sum_{j \in J} c_{ij} x_{ij} & \\ s.t. & \sum_{j \in J} x_{ij} =d_i & \forall i \in I \\ & \sum_{i \in I} x_{ij} \leq M_j y_j & \forall j \in J \\ & x_{ij} \leq d_i y_j & \forall i \in I; j \in J \\ & x_{ij} \geq 0 & \forall i \in I; j \in J \\ & y_j \in \{ 0,1 \} & \forall j \in J \end{array} \]
m = 2
n = 5
model = Model(name="Facility Location")
dummy_product = model.addProduct(name="dummy product")
supplier, customer = {}, {}
arc = {}
activity, mode ={}, {}
for i in range(m):
supplier[i] = model.addNode(name=f"Supplier({i})")
# model.addData(dtype="supply", product=dummy_product, node=supplier[i], under_penalty= 0, amount= 100 )
for j in range(n):
customer[j] = model.addNode(name=f"Customer({j})")
model.addData(dtype="demand", product=dummy_product, node=customer[j], amount= 20 )
for i in range(m):
for j in range(n):
arc[i,j] = model.addArc(name=f"arc({i},{j})", source=supplier[i], sink=customer[j])
activity[i,j] = model.addActivity(name=f"act({i},{j})", atype = "transport", product = dummy_product )
activity[i,j].addMode ( Mode(name=f"mode({i},{j})", variable_cost = random.randint(5,10) ) )
arc[i,j].addActivity( activity[i,j] )
#add resource and node activity
resource = {}
for i in range(m):
resource[i] = model.addResource(name=f"resource({i})", capacity = 100., fixed_cost=100.)
activity[i] = model.addActivity(name=f"act({i})",atype="make",product=dummy_product)
mode[i] = Mode(name=f"mode({i})", variable_cost = 1.)
mode[i].addResource(resource= resource[i], variable = 1.)
activity[i].addMode(mode[i])
supplier[i].addActivity(activity[i])
#pprint(model.model_dump_json(exclude=None))
gp_model = optimize(model)計画期間 = 0 1
Welcome to the CBC MILP Solver
Version: 2.10.3
Build Date: Dec 15 2019
command line - /Users/mikiokubo/Documents/dev/scmopt2/env/lib/python3.10/site-packages/pulp/solverdir/cbc/osx/64/cbc /var/folders/c0/1l3dlhys02nbkv4fyj9pkm680000gn/T/9f201fdedcfd49d19c61de181a863b07-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/c0/1l3dlhys02nbkv4fyj9pkm680000gn/T/9f201fdedcfd49d19c61de181a863b07-pulp.sol (default strategy 1)
At line 2 NAME MODEL
At line 3 ROWS
At line 26 COLUMNS
At line 139 RHS
At line 161 BOUNDS
At line 198 ENDATA
Problem MODEL has 21 rows, 36 columns and 60 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 820 - 0.00 seconds
Cgl0004I processed model has 7 rows, 22 columns (2 integer (2 of which binary)) and 32 elements
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 820
Cbc0038I Relaxing continuous gives 820
Cbc0038I Before mini branch and bound, 2 integers at bound fixed and 15 continuous
Cbc0038I Mini branch and bound did not improve solution (0.00 seconds)
Cbc0038I After 0.00 seconds - Feasibility pump exiting with objective of 820 - took 0.00 seconds
Cbc0012I Integer solution of 820 found by feasibility pump after 0 iterations and 0 nodes (0.00 seconds)
Cbc0001I Search completed - best objective 820, took 0 iterations and 0 nodes (0.00 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 820 to 820
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Result - Optimal solution found
Objective value: 820.00000000
Enumerated nodes: 0
Total iterations: 0
Time (CPU seconds): 0.00
Time (Wallclock seconds): 0.00
Option for printingOptions changed from normal to all
Total time (CPU seconds): 0.00 (Wallclock seconds): 0.00
Opt. Val= 820.0多品種流問題に枝上の固定費用 \(F: E \rightarrow \mathbf{R}_+\) をつけた問題を，多品種ネットワーク設計問題(multi-commodity network design problem)とよぶ． この問題は，\(NP\)-困難であり，しかも以下の通常の定式化だと小規模の問題例しか解けない．実務的には，パス型の定式化や列生成法を用いることが推奨される．
枝を使用するか否かを表す \(0\)-\(1\)変数を用いる．
目的関数： \[ \sum_{k \in K} \sum_{e \in E} c_e^k x_{e}^k + \sum_{e \in E} F_{ij} y_{ij} \]
フロー整合条件： \[ \sum_{j: ji \in E} x_{ji}^k - \sum_{j: ij \in E} x_{ij}^k = \left\{ \begin{array}{ll} -b_k & i=s_k \\ 0 & \forall i \in V \setminus \{s_k,t_k\} \\ b_k & i=t_k \end{array} \right. ,\forall k \in K \]
容量制約: \[ \sum_{k \in K} x_{ij}^k \leq u_{ij} y_{ij} \ \ \ \forall (i,j) \in E \]
非負制約: \[ x_{e}^k \geq 0 \ \ \ \forall e \in E, k \in K \]
整数制約： \[ y_{ij} \in \{ 0,1 \} \ \ \ \forall (i,j) \in E \]
random.seed(123)
m, n = 3, 3
cost_lb, cost_ub = 10, 10
cap_lb, cap_ub = 150, 150
demand_lb, demand_ub = 10, 30
G = nx.grid_2d_graph(m, n)
D = G.to_directed()
for (i, j) in D.edges():
D[i][j]["cost"] = random.randint(cost_lb, cost_ub)
D[i][j]["capacity"] = random.randint(cap_lb, cap_ub)
pos = {(i, j): (i, j) for (i, j) in G.nodes()}
b = {}
K = []
for i in D.nodes():
for j in D.nodes():
if i != j:
K.append((i, j))
b[i, j] = random.randint(demand_lb, demand_ub)
model = Model(name="Multicommodity Flow")
product = {}
node, arc = {}, {}
activity, mode ={}, {}
resource = {}
#node
for i in D.nodes():
node[i] = model.addNode(name=f"node({i})", location=pos[i])
#products
for (i,j) in b:
product[i,j] = model.addProduct(name=f"commodity{i},{j}")
model.addData(dtype="demand", product=product[i,j], node=node[j], amount=b[i,j])
model.addData(dtype="supply", product=product[i,j], node=node[i], amount=b[i,j])
#activity
for (i, j) in D.edges():
arc[i,j] = model.addArc(name=f"arc({i},{j})", source=node[i], sink=node[j])
resource[i,j] = model.addResource(name=f"resource({i},{j})", capacity=D[i][j]["capacity"], fixed_cost=100.)
for (k,l) in product:
activity[i,j,k,l] = model.addActivity(name=f"act({i},{j}.{k}.{l})", atype = "transport", product = product[k,l] )
mode[i,j,k,l] = Mode(name=f"mode({i},{j},{k},{l})", variable_cost = D[i][j]["cost"] )
mode[i,j,k,l].addResource(resource=resource[i,j], variable=1.)
activity[i,j,k,l].addMode ( mode[i,j,k,l] )
arc[i,j].addActivity( activity[i,j,k,l] )
model.update()
gp_model = optimize(model)計画期間 = 0 1
Welcome to the CBC MILP Solver
Version: 2.10.3
Build Date: Dec 15 2019
command line - /Users/mikiokubo/Documents/dev/scmopt2/env/lib/python3.10/site-packages/pulp/solverdir/cbc/osx/64/cbc /var/folders/c0/1l3dlhys02nbkv4fyj9pkm680000gn/T/77240d856e0c4406822031ba7c89b49c-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/c0/1l3dlhys02nbkv4fyj9pkm680000gn/T/77240d856e0c4406822031ba7c89b49c-pulp.sol (default strategy 1)
At line 2 NAME MODEL
At line 3 ROWS
At line 2405 COLUMNS
At line 16902 RHS
At line 19303 BOUNDS
At line 23072 ENDATA
Problem MODEL has 2400 rows, 3768 columns and 8952 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 30357.3 - 0.00 seconds
Cgl0004I processed model has 672 rows, 2040 columns (24 integer (24 of which binary)) and 5496 elements
Cbc0038I Initial state - 17 integers unsatisfied sum - 4.85333
Cbc0038I Pass 1: suminf. 0.00000 (0) obj. 3.54028e+08 iterations 321
Cbc0038I Solution found of 3.54028e+08
Cbc0038I Relaxing continuous gives 3.18026e+08
Cbc0038I Before mini branch and bound, 7 integers at bound fixed and 1891 continuous
Cbc0038I Full problem 672 rows 2040 columns, reduced to 26 rows 40 columns
Cbc0038I Mini branch and bound improved solution from 3.18026e+08 to 30860 (0.03 seconds)
Cbc0038I Round again with cutoff of 30809.7
Cbc0038I Pass 2: suminf. 1.13333 (3) obj. 30809.7 iterations 400
Cbc0038I Pass 3: suminf. 0.49733 (1) obj. 30809.7 iterations 241
Cbc0038I Pass 4: suminf. 0.42000 (1) obj. 30809.7 iterations 101
Cbc0038I Pass 5: suminf. 1.86667 (5) obj. 30809.7 iterations 201
Cbc0038I Pass 6: suminf. 1.45333 (4) obj. 30809.7 iterations 98
Cbc0038I Pass 7: suminf. 2.48000 (7) obj. 30809.7 iterations 153
Cbc0038I Pass 8: suminf. 1.80000 (5) obj. 30809.7 iterations 114
Cbc0038I Pass 9: suminf. 2.98000 (9) obj. 30809.7 iterations 211
Cbc0038I Pass 10: suminf. 1.89333 (5) obj. 30809.7 iterations 155
Cbc0038I Pass 11: suminf. 0.49733 (1) obj. 30809.7 iterations 211
Cbc0038I Pass 12: suminf. 0.42000 (1) obj. 30809.7 iterations 46
Cbc0038I Pass 13: suminf. 1.28667 (3) obj. 30809.7 iterations 195
Cbc0038I Pass 14: suminf. 1.28667 (3) obj. 30809.7 iterations 0
Cbc0038I Pass 15: suminf. 1.62000 (5) obj. 30809.7 iterations 159
Cbc0038I Pass 16: suminf. 1.61333 (4) obj. 30809.7 iterations 35
Cbc0038I Pass 17: suminf. 2.28000 (7) obj. 30809.7 iterations 107
Cbc0038I Pass 18: suminf. 1.18667 (3) obj. 30809.7 iterations 141
Cbc0038I Pass 19: suminf. 0.49733 (1) obj. 30809.7 iterations 118
Cbc0038I Pass 20: suminf. 0.42000 (1) obj. 30809.7 iterations 64
Cbc0038I Pass 21: suminf. 2.25333 (6) obj. 30809.7 iterations 286
Cbc0038I Pass 22: suminf. 1.72666 (5) obj. 30809.7 iterations 103
Cbc0038I Pass 23: suminf. 2.78666 (8) obj. 30809.7 iterations 141
Cbc0038I Pass 24: suminf. 2.38666 (6) obj. 30809.7 iterations 136
Cbc0038I Pass 25: suminf. 1.42000 (4) obj. 30809.7 iterations 80
Cbc0038I Pass 26: suminf. 0.49733 (1) obj. 30809.7 iterations 113
Cbc0038I Pass 27: suminf. 0.42000 (1) obj. 30809.7 iterations 44
Cbc0038I Pass 28: suminf. 2.23333 (6) obj. 30809.7 iterations 292
Cbc0038I Pass 29: suminf. 1.70000 (4) obj. 30809.7 iterations 74
Cbc0038I Pass 30: suminf. 2.31333 (6) obj. 30809.7 iterations 54
Cbc0038I Pass 31: suminf. 2.14667 (5) obj. 30809.7 iterations 3
Cbc0038I No solution found this major pass
Cbc0038I Before mini branch and bound, 3 integers at bound fixed and 1775 continuous
Cbc0038I Full problem 672 rows 2040 columns, reduced to 71 rows 106 columns
Cbc0038I Mini branch and bound did not improve solution (0.14 seconds)
Cbc0038I After 0.14 seconds - Feasibility pump exiting with objective of 30860 - took 0.12 seconds
Cbc0012I Integer solution of 30860 found by feasibility pump after 0 iterations and 0 nodes (0.15 seconds)
Cbc0038I Full problem 672 rows 2040 columns, reduced to 672 rows 2033 columns - 14 fixed gives 672, 2019 - still too large
Cbc0031I 25 added rows had average density of 11.64
Cbc0013I At root node, 42 cuts changed objective from 30357.333 to 30860 in 6 passes
Cbc0014I Cut generator 0 (Probing) - 76 row cuts average 2.0 elements, 0 column cuts (0 active) in 0.002 seconds - new frequency is 1
Cbc0014I Cut generator 1 (Gomory) - 37 row cuts average 317.1 elements, 0 column cuts (0 active) in 0.002 seconds - new frequency is 1
Cbc0014I Cut generator 2 (Knapsack) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.003 seconds - new frequency is -100
Cbc0014I Cut generator 3 (Clique) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100
Cbc0014I Cut generator 4 (MixedIntegerRounding2) - 55 row cuts average 9.4 elements, 0 column cuts (0 active) in 0.001 seconds - new frequency is 1
Cbc0014I Cut generator 5 (FlowCover) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.002 seconds - new frequency is -100
Cbc0014I Cut generator 6 (TwoMirCuts) - 64 row cuts average 32.4 elements, 0 column cuts (0 active) in 0.002 seconds - new frequency is 1
Cbc0001I Search completed - best objective 30860, took 343 iterations and 0 nodes (0.18 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 30357.3 to 30860
Probing was tried 6 times and created 76 cuts of which 0 were active after adding rounds of cuts (0.002 seconds)
Gomory was tried 6 times and created 37 cuts of which 0 were active after adding rounds of cuts (0.002 seconds)
Knapsack was tried 6 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.003 seconds)
Clique was tried 6 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 6 times and created 55 cuts of which 0 were active after adding rounds of cuts (0.001 seconds)
FlowCover was tried 6 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.002 seconds)
TwoMirCuts was tried 6 times and created 64 cuts of which 0 were active after adding rounds of cuts (0.002 seconds)
ZeroHalf was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Result - Optimal solution found
Objective value: 30860.00000000
Enumerated nodes: 0
Total iterations: 343
Time (CPU seconds): 0.19
Time (Wallclock seconds): 0.19
Option for printingOptions changed from normal to all
Total time (CPU seconds): 0.19 (Wallclock seconds): 0.20
Opt. Val= 30860.0多期間複数ラインのロジスティクス・ネットワーク設計モデルの例題である。某飲料メーカーの事例をもとにしている。
JSON VIsualizer https://omute.net/ を使う。
cost (Node,Arc,Mode)にfixed,variable,piecewiseを定義
変数
制約
\[ I_{j,t-1}^p + \sum_{i,m} X_{ijmt}^p = \sum_{j,m} d_{jt}^p + I_{i,t}^p \ \ \ \forall j,p,t \]
\[ \sum_{m,p} varriable_{ijm} X_{imt}^p + \sum_{m,p} setup_{ijm} Y_{imt}^p \leq CAP_{ijrt} \ \ \ \forall i,j,r,t \]
\[ x_{imt}^p \leq M y_{imt}^p \ \ \ \forall i,m,t,p \]
\[ X_{ijmt}^p \leq M Y_{ijmt}^p \ \ \ \forall i,j,m,t,p \]
TODO: 一般化
Welcome to the CBC MILP Solver
Version: 2.10.3
Build Date: Dec 15 2019
command line - /Users/mikiokubo/Documents/dev/scmopt2/env/lib/python3.10/site-packages/pulp/solverdir/cbc/osx/64/cbc /var/folders/c0/1l3dlhys02nbkv4fyj9pkm680000gn/T/4421000aad3442f3b08b6b5f90380a78-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/c0/1l3dlhys02nbkv4fyj9pkm680000gn/T/4421000aad3442f3b08b6b5f90380a78-pulp.sol (default strategy 1)
At line 2 NAME MODEL
At line 3 ROWS
At line 71 COLUMNS
At line 436 RHS
At line 503 BOUNDS
At line 588 ENDATA
Problem MODEL has 66 rows, 84 columns and 232 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 240.06 - 0.00 seconds
Cgl0003I 0 fixed, 0 tightened bounds, 24 strengthened rows, 0 substitutions
Cgl0004I processed model has 66 rows, 84 columns (24 integer (24 of which binary)) and 232 elements
Cbc0038I Initial state - 9 integers unsatisfied sum - 0.618557
Cbc0038I Pass 1: suminf. 0.61856 (3) obj. 304.856 iterations 18
Cbc0038I Solution found of 543
Cbc0038I Relaxing continuous gives 543
Cbc0038I Before mini branch and bound, 12 integers at bound fixed and 26 continuous
Cbc0038I Full problem 66 rows 84 columns, reduced to 19 rows 28 columns
Cbc0038I Mini branch and bound did not improve solution (0.01 seconds)
Cbc0038I Round again with cutoff of 518.886
Cbc0038I Pass 2: suminf. 0.61856 (3) obj. 304.856 iterations 0
Cbc0038I Pass 3: suminf. 0.24114 (1) obj. 518.886 iterations 5
Cbc0038I Pass 4: suminf. 0.00000 (0) obj. 446 iterations 5
Cbc0038I Solution found of 446
Cbc0038I Relaxing continuous gives 446
Cbc0038I Before mini branch and bound, 12 integers at bound fixed and 25 continuous
Cbc0038I Full problem 66 rows 84 columns, reduced to 17 rows 25 columns
Cbc0038I Mini branch and bound did not improve solution (0.01 seconds)
Cbc0038I Round again with cutoff of 417.171
Cbc0038I Pass 5: suminf. 0.61856 (3) obj. 304.856 iterations 0
Cbc0038I Pass 6: suminf. 0.67066 (2) obj. 417.171 iterations 4
Cbc0038I Pass 7: suminf. 0.28829 (1) obj. 417.171 iterations 3
Cbc0038I Pass 8: suminf. 0.30928 (1) obj. 376.928 iterations 2
Cbc0038I Pass 9: suminf. 0.68566 (3) obj. 417.171 iterations 8
Cbc0038I Pass 10: suminf. 0.68566 (3) obj. 417.171 iterations 3
Cbc0038I Pass 11: suminf. 0.69566 (2) obj. 417.171 iterations 4
Cbc0038I Pass 12: suminf. 0.49447 (2) obj. 417.171 iterations 3
Cbc0038I Pass 13: suminf. 0.28829 (1) obj. 417.171 iterations 10
Cbc0038I Pass 14: suminf. 0.30928 (1) obj. 376.928 iterations 1
Cbc0038I Pass 15: suminf. 0.76671 (4) obj. 417.171 iterations 19
Cbc0038I Pass 16: suminf. 0.30928 (1) obj. 373.928 iterations 9
Cbc0038I Pass 17: suminf. 0.25829 (1) obj. 417.171 iterations 2
Cbc0038I Pass 18: suminf. 0.76671 (5) obj. 417.171 iterations 19
Cbc0038I Pass 19: suminf. 0.51546 (3) obj. 393.546 iterations 6
Cbc0038I Pass 20: suminf. 0.75171 (3) obj. 417.171 iterations 8
Cbc0038I Pass 21: suminf. 0.77171 (5) obj. 417.171 iterations 10
Cbc0038I Pass 22: suminf. 0.76171 (4) obj. 417.171 iterations 15
Cbc0038I Pass 23: suminf. 0.65566 (3) obj. 417.171 iterations 4
Cbc0038I Pass 24: suminf. 0.65566 (3) obj. 417.171 iterations 5
Cbc0038I Pass 25: suminf. 0.66566 (2) obj. 417.171 iterations 6
Cbc0038I Pass 26: suminf. 0.46447 (2) obj. 417.171 iterations 3
Cbc0038I Pass 27: suminf. 0.25829 (1) obj. 417.171 iterations 8
Cbc0038I Pass 28: suminf. 0.30928 (1) obj. 373.928 iterations 1
Cbc0038I Pass 29: suminf. 0.76671 (4) obj. 417.171 iterations 12
Cbc0038I Pass 30: suminf. 0.30928 (1) obj. 373.928 iterations 13
Cbc0038I Pass 31: suminf. 0.76671 (4) obj. 417.171 iterations 10
Cbc0038I Pass 32: suminf. 0.51546 (3) obj. 393.046 iterations 12
Cbc0038I Pass 33: suminf. 0.75671 (3) obj. 417.171 iterations 4
Cbc0038I Pass 34: suminf. 0.75671 (4) obj. 417.171 iterations 20
Cbc0038I Rounding solution of 443 is better than previous of 446
Cbc0038I Before mini branch and bound, 5 integers at bound fixed and 9 continuous
Cbc0038I Full problem 66 rows 84 columns, reduced to 45 rows 69 columns
Cbc0038I Mini branch and bound did not improve solution (0.01 seconds)
Cbc0038I Round again with cutoff of 382.576
Cbc0038I Pass 34: suminf. 0.61856 (3) obj. 304.856 iterations 0
Cbc0038I Pass 35: suminf. 0.63516 (3) obj. 382.576 iterations 3
Cbc0038I Pass 36: suminf. 0.30928 (1) obj. 376.928 iterations 7
Cbc0038I Pass 37: suminf. 0.36576 (1) obj. 382.576 iterations 4
Cbc0038I Pass 38: suminf. 0.41398 (5) obj. 382.576 iterations 11
Cbc0038I Pass 39: suminf. 0.30928 (1) obj. 376.928 iterations 9
Cbc0038I Pass 40: suminf. 0.36576 (1) obj. 382.576 iterations 2
Cbc0038I Pass 41: suminf. 0.41898 (4) obj. 382.576 iterations 6
Cbc0038I Pass 42: suminf. 0.41237 (4) obj. 382.576 iterations 5
Cbc0038I Pass 43: suminf. 0.42076 (4) obj. 382.576 iterations 2
Cbc0038I Pass 44: suminf. 0.42076 (4) obj. 382.576 iterations 8
Cbc0038I Pass 45: suminf. 0.61016 (5) obj. 382.576 iterations 10
Cbc0038I Pass 46: suminf. 0.42076 (4) obj. 382.576 iterations 7
Cbc0038I Pass 47: suminf. 0.41237 (3) obj. 382.576 iterations 7
Cbc0038I Pass 48: suminf. 0.41576 (3) obj. 382.576 iterations 1
Cbc0038I Pass 49: suminf. 0.61016 (5) obj. 382.576 iterations 11
Cbc0038I Pass 50: suminf. 0.30928 (1) obj. 373.928 iterations 7
Cbc0038I Pass 51: suminf. 0.39576 (1) obj. 382.576 iterations 2
Cbc0038I Pass 52: suminf. 0.61516 (4) obj. 382.576 iterations 4
Cbc0038I Pass 53: suminf. 0.60516 (6) obj. 382.576 iterations 8
Cbc0038I Pass 54: suminf. 0.62016 (5) obj. 382.576 iterations 4
Cbc0038I Pass 55: suminf. 0.30928 (1) obj. 376.928 iterations 11
Cbc0038I Pass 56: suminf. 0.36576 (1) obj. 382.576 iterations 3
Cbc0038I Pass 57: suminf. 0.41576 (3) obj. 382.576 iterations 13
Cbc0038I Pass 58: suminf. 0.41237 (3) obj. 382.576 iterations 2
Cbc0038I Pass 59: suminf. 0.41576 (3) obj. 382.576 iterations 1
Cbc0038I Pass 60: suminf. 0.61016 (5) obj. 382.576 iterations 7
Cbc0038I Pass 61: suminf. 0.42076 (4) obj. 382.576 iterations 11
Cbc0038I Pass 62: suminf. 0.41898 (4) obj. 382.576 iterations 15
Cbc0038I Pass 63: suminf. 0.30928 (1) obj. 373.928 iterations 6
Cbc0038I No solution found this major pass
Cbc0038I Before mini branch and bound, 2 integers at bound fixed and 8 continuous
Cbc0038I Full problem 66 rows 84 columns, reduced to 41 rows 64 columns
Cbc0038I Mini branch and bound did not improve solution (0.02 seconds)
Cbc0038I After 0.02 seconds - Feasibility pump exiting with objective of 443 - took 0.01 seconds
Cbc0012I Integer solution of 443 found by feasibility pump after 0 iterations and 0 nodes (0.02 seconds)
Cbc0038I Full problem 66 rows 84 columns, reduced to 43 rows 52 columns
Cbc0031I 10 added rows had average density of 27.8
Cbc0013I At root node, 10 cuts changed objective from 301.85567 to 442.78042 in 19 passes
Cbc0014I Cut generator 0 (Probing) - 1 row cuts average 22.0 elements, 0 column cuts (0 active) in 0.001 seconds - new frequency is -100
Cbc0014I Cut generator 1 (Gomory) - 77 row cuts average 29.0 elements, 0 column cuts (0 active) in 0.001 seconds - new frequency is 1
Cbc0014I Cut generator 2 (Knapsack) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.002 seconds - new frequency is -100
Cbc0014I Cut generator 3 (Clique) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100
Cbc0014I Cut generator 4 (MixedIntegerRounding2) - 20 row cuts average 22.4 elements, 0 column cuts (0 active) in 0.001 seconds - new frequency is 1
Cbc0014I Cut generator 5 (FlowCover) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.001 seconds - new frequency is -100
Cbc0014I Cut generator 6 (TwoMirCuts) - 42 row cuts average 20.3 elements, 0 column cuts (0 active) in 0.001 seconds - new frequency is 1
Cbc0010I After 0 nodes, 1 on tree, 443 best solution, best possible 442.78042 (0.03 seconds)
Cbc0001I Search completed - best objective 443, took 322 iterations and 2 nodes (0.03 seconds)
Cbc0032I Strong branching done 18 times (116 iterations), fathomed 2 nodes and fixed 2 variables
Cbc0035I Maximum depth 0, 8 variables fixed on reduced cost
Cuts at root node changed objective from 301.856 to 442.78
Probing was tried 19 times and created 1 cuts of which 0 were active after adding rounds of cuts (0.001 seconds)
Gomory was tried 21 times and created 77 cuts of which 0 were active after adding rounds of cuts (0.001 seconds)
Knapsack was tried 19 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.002 seconds)
Clique was tried 19 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 21 times and created 20 cuts of which 0 were active after adding rounds of cuts (0.001 seconds)
FlowCover was tried 19 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.001 seconds)
TwoMirCuts was tried 21 times and created 42 cuts of which 0 were active after adding rounds of cuts (0.001 seconds)
ZeroHalf was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Result - Optimal solution found
Objective value: 443.00000000
Enumerated nodes: 2
Total iterations: 322
Time (CPU seconds): 0.03
Time (Wallclock seconds): 0.03
Option for printingOptions changed from normal to all
Total time (CPU seconds): 0.03 (Wallclock seconds): 0.03
print("Obj. Val=",gp_model.ObjVal)
for (i,m,p,t) in x:
if x[i,m,p,t].X > 0:
print(i,m,t,p,x[i,m,p,t].X)Obj. Val= 443.00000160000013
plant(1) mode(1,1) 0 prod(0) 30.0
plant(1) mode(1,1) 0 prod(1) 30.0plant(1) customer(0) trans_mode(1,0,0) 0 prod(0) 11.666667
plant(1) customer(0) trans_mode(1,0,0) 1 prod(0) 1.6666667
plant(1) customer(0) trans_mode(1,0,0) 2 prod(0) 1.6666667
plant(1) customer(0) trans_mode(1,0,1) 0 prod(1) 11.666667
plant(1) customer(0) trans_mode(1,0,1) 1 prod(1) 1.6666667
plant(1) customer(0) trans_mode(1,0,1) 2 prod(1) 1.6666667
plant(1) customer(1) trans_mode(1,1,0) 0 prod(0) 11.666667
plant(1) customer(1) trans_mode(1,1,0) 1 prod(0) 1.6666667
plant(1) customer(1) trans_mode(1,1,0) 2 prod(0) 1.6666667
plant(1) customer(1) trans_mode(1,1,1) 0 prod(1) 11.666667
plant(1) customer(1) trans_mode(1,1,1) 1 prod(1) 1.6666667
plant(1) customer(1) trans_mode(1,1,1) 2 prod(1) 1.6666667ロットサイズ決定モデルの例題である。某化学メーカーの事例をもとにしている。
集合:
パラメータ:
変数:
上の記号を用いると、多段階ロットサイズ決定モデルは，以下のように定式化できる．
\[ \begin{array}{ l l l } minimize & \sum_{t=1}^T \sum_{p \in P} \left( f_t^p y_t^p + c_t^p x_t^p + h_t^p I_t^p \right) & \\ s.t. & \ \ I_{t-1}^p +x_t^p = d_t^p+ \sum_{q \in Parent_p} \phi_{pq} x_t^q +I_t^p & \forall p \in P, t=1,\ldots,T \\ & \sum_{p \in P_k} x_t^p +\sum_{p \in P_k} g_t^p y_t^p \leq M_t^k & \forall k \in K, t=1,\ldots,T \\ & x_t^p \leq UB_t^p y_t^p & \forall p \in P, t=1,\ldots,T \\ & I_0^p =0 & \forall p \in P \\ & x_t^p,I_t^p \geq 0 & \forall p \in P, t=1,\ldots,T \\ & y_t \in \{0,1\} & \forall t=1,\ldots,T \end{array} \]
上の定式化で，最初の制約式は，各期および各品目に対する在庫の保存式を表す． より具体的には，品目 \(p\) の期 \(t-1\) からの在庫量 \(I_{t-1}^p\) と生産量 \(x_t^p\) を加えたものが， 期 \(t\) における需要量 \(d_t^p\)，次期への在庫量 \(I_t^p\)， および他の品目を生産するときに必要な量 \(\sum_{q \in Parent_p} \phi_{pq} x_t^q\) の和に等しいことを表す．
2番目の制約は， 各期の生産時間の上限制約を表す． 定式化ではすべての品目の生産時間は， \(1\) 単位あたり\(1\) 時間になるようにスケーリングしてあると仮定していたが，実際問題のモデル化の際には， 品目 p を \(1\) 単位生産されるときに，資源 \(r\) を使用する時間を用いた方が汎用性がある．
3番目の式は，段取り替えをしない期は生産できないことを表す．
prod_df = pd.read_csv(folder + "lotprod.csv",index_col=0)
prod_df.set_index("name", inplace=True)
production_df = pd.read_csv(folder + "production.csv",index_col=0)
bom_df = pd.read_csv(folder + "bomodel.csv", index_col =0)
resource_df = pd.read_csv(folder + "resource.csv", index_col=0)
plnt_demand_df = pd.read_csv(folder+"plnt-demand.csv", index_col=0)
demand_df = pd.pivot_table(plnt_demand_df, index= "prod", columns ="period", values="demand", aggfunc="sum")
raw_materials = set(bom_df["child"])
final_products = set(bom_df["parent"])
items = raw_materials | final_products| inv_cost | safety_inventory | initial_inventory | target_inventory | |
|---|---|---|---|---|
| name | ||||
| A | 0.005753 | 1236.729760 | 4972.0 | 8708.408003 |
| B | 0.005753 | 1255.397379 | 4945.0 | 8636.582305 |
| C | 0.005753 | 499.358809 | 3149.0 | 5798.853925 |
| D | 0.005753 | 33.835059 | 908.0 | 1782.369511 |
| E | 0.005753 | 5653.955112 | 14448.0 | 23243.815827 |
| F | 0.005753 | 847.043213 | 4083.0 | 7320.560359 |
| G | 0.005753 | 1197.061070 | 5268.0 | 9339.456295 |
| H | 0.005753 | 2872.479877 | 7541.0 | 12211.145021 |
| I | 0.005753 | 168.008571 | 1703.0 | 3239.099583 |
| J | 0.005753 | 3788.359936 | 9955.0 | 16123.176759 |
| ABCDE | 0.002000 | 8679.276120 | 28422.0 | 48170.029571 |
| FGHIJ | 0.002000 | 8872.952667 | 28550.0 | 48233.438016 |
| name | ProdTime | SetupTime | ProdCost | SetupCost | |
|---|---|---|---|---|---|
| 0 | A | 1 | 4815 | 263 | 13922 |
| 1 | B | 1 | 3820 | 237 | 12045 |
| 2 | C | 1 | 6687 | 235 | 17853 |
| 3 | D | 1 | 4892 | 147 | 13550 |
| 4 | E | 1 | 7121 | 190 | 11467 |
| 5 | F | 1 | 6241 | 149 | 11939 |
| 6 | G | 1 | 5868 | 154 | 14069 |
| 7 | H | 1 | 4913 | 176 | 19382 |
| 8 | I | 1 | 3682 | 257 | 10532 |
| 9 | J | 1 | 3696 | 187 | 19639 |
| 10 | ABCDE | 1 | 6002 | 202 | 16403 |
| 11 | FGHIJ | 1 | 6644 | 151 | 12775 |
| child | parent | units | |
|---|---|---|---|
| 0 | ABCDE | A | 1 |
| 1 | ABCDE | B | 1 |
| 2 | ABCDE | C | 1 |
| 3 | ABCDE | D | 1 |
| 4 | ABCDE | E | 1 |
| 5 | FGHIJ | F | 1 |
| 6 | FGHIJ | G | 1 |
| 7 | FGHIJ | H | 1 |
| 8 | FGHIJ | I | 1 |
| 9 | FGHIJ | J | 1 |
| prod | period | demand | |
|---|---|---|---|
| 0 | A | 0 | 12771 |
| 1 | A | 1 | 11228 |
| 2 | A | 2 | 10151 |
| 3 | A | 3 | 11464 |
| 4 | A | 4 | 12076 |
| ... | ... | ... | ... |
| 235 | J | 19 | 46408 |
| 236 | J | 20 | 52901 |
| 237 | J | 21 | 53650 |
| 238 | J | 22 | 49522 |
| 239 | J | 23 | 47466 |
240 rows × 3 columns
| period | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| prod | |||||||||||||||||||||
| A | 12771 | 11228 | 10151 | 11464 | 12076 | 9687 | 11828 | 11361 | 12026 | 11132 | ... | 14854 | 17207 | 17868 | 15403 | 16722 | 15307 | 17087 | 14123 | 16964 | 15739 |
| B | 11352 | 11820 | 11174 | 11515 | 10768 | 10059 | 10928 | 11325 | 11565 | 10215 | ... | 15937 | 15124 | 15590 | 16225 | 15623 | 15478 | 16595 | 16567 | 15249 | 16403 |
| C | 4212 | 4448 | 4864 | 3999 | 4570 | 4189 | 4542 | 4784 | 4019 | 4522 | ... | 5623 | 5845 | 6421 | 6910 | 6187 | 6035 | 6511 | 5725 | 6559 | 6701 |
| D | 497 | 488 | 460 | 525 | 516 | 508 | 532 | 482 | 510 | 528 | ... | 748 | 736 | 754 | 756 | 761 | 717 | 730 | 761 | 715 | 731 |
| E | 50189 | 47903 | 52992 | 52302 | 47775 | 48568 | 52320 | 46079 | 49140 | 51636 | ... | 69792 | 67832 | 66657 | 75566 | 65866 | 70182 | 65074 | 70326 | 75577 | 72824 |
| F | 7715 | 7850 | 8021 | 7018 | 7618 | 7797 | 6750 | 7872 | 7189 | 7054 | ... | 10281 | 10440 | 11190 | 10920 | 9685 | 10489 | 10707 | 10287 | 11073 | 10841 |
| G | 9666 | 10178 | 9787 | 10356 | 9343 | 10045 | 9861 | 10378 | 10249 | 9379 | ... | 14998 | 15079 | 13523 | 12798 | 14650 | 13749 | 13587 | 13978 | 13924 | 15514 |
| H | 22833 | 24903 | 27112 | 25403 | 27764 | 25263 | 26469 | 22517 | 28560 | 25838 | ... | 33085 | 35007 | 33663 | 32769 | 37956 | 38661 | 35197 | 36849 | 36031 | 37895 |
| I | 1589 | 1469 | 1570 | 1677 | 1575 | 1610 | 1553 | 1756 | 1621 | 1580 | ... | 2314 | 2436 | 2363 | 2343 | 2611 | 2175 | 2214 | 2145 | 1897 | 2276 |
| J | 32025 | 32982 | 33393 | 36982 | 37142 | 30692 | 34577 | 33559 | 34170 | 32205 | ... | 44326 | 52677 | 48896 | 49668 | 49541 | 46408 | 52901 | 53650 | 49522 | 47466 |
10 rows × 24 columns
model = Model(name="lot sizing model")
model.interest_rate = 0.15
T = len(demand_df.columns)
periods, products, resources, arcs = {}, {}, {}, {}
activities, modes = {}, {}
#period
for t in range(T):
periods[t] = model.addPeriod(name=t)
#products
for row in prod_df.itertuples():
p = str(row.Index)
products[p] = model.addProduct(name=f"prod({p})", value=row.inv_cost/model.interest_rate)
#demand
demand_matrix = demand_df.values
for i,p in enumerate(products):
if i>=10:
break #部品の需要はなし
for t in periods:
model.addData(dtype="demand",product=products[p],period=periods[t],amount= float(demand_matrix[i,t]))
#resource 期ごとに同じ容量を仮定（２つの資源があり，Res0は原料生産，Res1は最終製品生産に用いられる）
for row in resource_df.itertuples():
resources[row.name] = model.addResource(name=row.name, capacity=row.capacity)
#activity and mode
for row in production_df.itertuples():
p = row.name
modes[p] = model.addMode(name=p, fixed_cost=row.SetupCost, variable_cost=row.ProdCost)
if p in raw_materials:
modes[p].addResource(resources["Res0"], fixed=row.SetupTime, variable = row.ProdTime)
elif p in final_products:
modes[p].addResource(resources["Res1"], fixed=row.SetupTime, variable = row.ProdTime)
activities[p] = model.addActivity(name=p, atype="make", product=products[p])
activities[p].addMode( modes[p] )
#bom
for row in bom_df.itertuples():
child = row.child
p = row.parent
modes[p].addComponent(products[child])リード時間 \(0\) を仮定
変数
制約
\[ I_{t-1}^p + \sum_{m} x_{mt}^p = d_{t}^p + \sum_{q} \phi_{pqm} x_{mt}^q + I_{i,t}^p \ \ \ \forall t,p \]
\[ \sum_{m,p} varriable_{m} x_{mt}^p + \sum_{m,p} setup_{m} y_{mt}^p \leq CAP_{rt} \ \ \ \forall r,t \]
\[ x_{mt}^p \leq M y_{mt}^p \ \ \ \forall m,t,p \]
print("Obj. Val.=",gp_model.ObjVal)
for i in range(5):
print(cost[i].X)
# for (t,m,p) in x:
# if x[t,m,p].X > 0:
# print(t,m,p, x[t,m,p].X)Obj. Val.= 21943968879553.82
2194302.0
2185400.0
951069700.0
4155.822
None0 prod(ABCDE) 18377.5
1 prod(ABCDE) 39889.0
2 prod(ABCDE) 57646.5
3 prod(ABCDE) 75240.0
5 prod(ABCDE) 24387.5
6 prod(ABCDE) 41636.0
7 prod(ABCDE) 65003.5
8 prod(ABCDE) 85142.0
9 prod(ABCDE) 104507.5
10 prod(ABCDE) 122496.0
11 prod(ABCDE) 140690.5
12 prod(ABCDE) 129483.0
13 prod(ABCDE) 110703.5
14 prod(ABCDE) 101148.0
15 prod(ABCDE) 91802.5
16 prod(ABCDE) 81911.0
17 prod(ABCDE) 64449.5
18 prod(ABCDE) 56689.0
19 prod(ABCDE) 46368.5
20 prod(ABCDE) 37770.0
21 prod(ABCDE) 27666.5
22 prod(ABCDE) 10001.0
4 prod(FGHIJ) 89138.5
5 prod(FGHIJ) 81341.5
6 prod(FGHIJ) 74591.5
7 prod(FGHIJ) 66719.5
8 prod(FGHIJ) 59530.5
9 prod(FGHIJ) 52476.5
10 prod(FGHIJ) 44017.5
11 prod(FGHIJ) 36125.5
12 prod(FGHIJ) 24921.5
13 prod(FGHIJ) 13160.5
14 prod(FGHIJ) 2879.5簡単なシフト最適化モデル
TODO: 複数日の問題への拡張
変数
目的関数
\[ minimize \ \ \ \sum_{a,m} cost_{am} x_{am} \]
制約
\[ \sum_{m} x_{am} \leq 1 \ \ \ \forall a \]
\[ \sum_{a,m | t \in period[a,m] } x_{am} \geq demand_{t} \ \ \ \forall t \]
日をまたぐ制約をConstraintクラスで記述
0 period(0,0) 1 1
1 period(0,1) 2 2
2 period(0,2) 3 3
3 period(0,3) 4 4
4 period(0,4) 5 5
5 period(0,5) 4 4
6 period(0,6) 3 3
7 period(0,7) 2 2
8 period(0,8) 2 2
9 period(0,9) 1 1
10 period(1,0) 1 1
11 period(1,1) 2 2
12 period(1,2) 3 3
13 period(1,3) 4 4
14 period(1,4) 5 5
15 period(1,5) 5 5
16 period(1,6) 3 3
17 period(1,7) 2 2
18 period(1,8) 2 2
19 period(1,9) 1 1
20 period(2,0) 2 2
21 period(2,1) 2 2
22 period(2,2) 3 3
23 period(2,3) 4 4
24 period(2,4) 5 5
25 period(2,5) 5 5
26 period(2,6) 3 3
27 period(2,7) 2 2
28 period(2,8) 1 1
29 period(2,9) 1 1枝上に活動を定義した一般化ロジスティクス・ネットワーク設計モデル
可視化の際に枝の太さでフロー量を表現できる．
注：区分的線形関数をSOSで解く
集合：
\(N\): 点の集合．原料供給地点，工場，倉庫の配置可能地点，顧客群の集合などのすべての地点を総称して点とよぶ．
\(A\): 枝の集合． 少なくとも1つの製品が移動する可能性のある点の対を枝とよぶ．
\(Prod\): 製品の集合． 製品は，ロジスティクス・ネットワーク内を流れる「もの」の総称である．
以下に定義する \(Child_p\)， \(Parent_p\) は，製品の集合の部分集合である．
\(Child_p\): 部品展開表における製品 \(p\) の子製品の集合．言い換えれば，製品 \(p\) を製造するために必要な製品の集合．
\(Parent_p\): 部品展開表における製品 \(p\) の親製品の集合．言い換えれば，製品 \(p\) を分解することによって生成される製品の集合．
\(Res\): 資源の集合． 製品を処理（移動，組み立て，分解）するための資源の総称． 基本モデルでは枝上で定義される． たとえば，工場を表す枝における生産ライン（もしくは機械）や 輸送を表す枝における輸送機器（トラック，船，鉄道，飛行機など）が資源の代表的な要素となる．
\(NodeProd\): 需要もしくは供給が発生する点と製品の \(2\)つ組の集合．
\(ArcRes\): 枝と資源の可能な \(2\)つ組の集合． 枝 \(a \in A\) 上で資源 \(r \in Res\) が利用可能なとき，\((a,r)\) の組が 集合 \(ArcRes\) に含まれるものとする．
\(ArcResProd\): 枝と資源と製品の可能な \(3\)つ組の集合． 枝 \(a \in A\) 上の資源 \(r \in Res\) で製品 \(p \in Prod\) の処理が利用可能なとき， \((a,r,p)\) の組が 集合 \(ArcResProd\) に含まれるものとする．
以下に定義する \(Assemble\)，\(Disassemble\) および \(Transit\) は \(ArcResProd\) の部分集合である．
\(Assemble\): 組み立てを表す枝と資源と製品の可能な \(3\) つ組の集合． 枝 \(a \in A\) 上の資源 \(r \in Res\) で製品 \(p \in Prod\) の組み立て処理が利用可能なとき，\((a,r,p)\) の組が 集合 \(Assemble\) に含まれるものとする．ここで製品 \(p\) の組み立て処理とは，子製品の集合 \(Child_p\) を用いて \(p\) を 製造することを指す．
\(Disassemble\): 分解を表す枝と資源と製品の可能な \(3\) つ組の集合． 枝 \(a \in A\) 上の資源 \(r \in Res\) で製品 \(p \in Prod\) の分解処理が利用可能なとき，\((a,r,p)\) の組が 集合 \(Disassemble\) に含まれるものとする．ここで製品 \(p\) の分解処理とは，\(p\) から親製品の集合 \(Parent_p\) を 生成することを指す．
\(Transit\): 移動を表す枝と資源と製品の可能な \(3\) つ組の集合． 枝 \(a \in A\) 上の資源 \(r \in Res\) で製品 \(p \in Prod\) が形態を変えずに流れることが可能なとき， \((a,r,p)\) の組は集合 \(Transit\) に含まれるものとする．
\(ResProd\): 資源と製品の可能な \(2\) つ組の集合． 集合 \(ArcResProd\) から生成される．
\(ArcProd\): 枝と製品の可能な \(2\) つ組の集合． 集合 \(ArcResProd\) から生成される．
入力データ：
\(D_i^p\): 点 \(i\) における製品 \(p\) の需要量（\(p\)-units \(/\) 単位期間）； 負の需要は供給量を表す．ここで，\(p\)-unit とは，製品 \(p\) の \(1\) 単位を表す．
\(DPENALTY_{ip}^{+}\): 点 \(i\) における製品 \(p\) の \(1\) 単位あたりの需要超過（供給余裕）ペナルティ （円 \(/\) 単位期間・\(p\)-unit）；通常は小さな値
\(DPENALTY_{ip}^{-}\): 点 \(i\) における製品 \(p\) の \(1\) 単位あたりの需要不足（供給超過）ペナルティ （円 \(/\) 単位期間・\(p\)-unit）； 通常は大きな値
\(AFC_{ij}\): 枝 \((i,j)\) を使用するときに発生する固定費用（円 \(/\) 単位期間）
\(ARFC_{ijr}:\) 枝 \((i,j)\) 上で資源 \(r\) を使用するときに発生する固定費用（円 \(/\) 単位期間）
\(ARPVC_{ijr}^p\): 枝 \((i,j)\) 上で資源 \(r\) を利用して製品 \(p\) を \(1\) 単位処理するごとに発生する変動費用 （円 \(/\) 単位期間・\(p\)-unit）
\(\phi_{pq}\) : \(q \in Parent_p\) のとき， 品目 \(q\) を \(1\) 単位生成するのに必要な品目 \(p\) の数 （\(p\)-units）； ここで， \(p\)-unitsとは，品目 \(q\) の \(1\)単位と混同しないために導入された単位であり， 品目 \(p\) の \(1\)単位を表す．\(\phi_{pq}\) は，部品展開表を有向グラフ表現したときには，枝の重みを表す． この値から以下の\(U_{p q}\)と\(\bar{U}_{p q}\)を計算する．
\(U_{p q}\): 製品 \(p\) の \(1\) 単位を組み立て処理するために必要な製品 \(q \in Child_p\) の量（\(q\)-units）
\(\bar{U}_{p q}\): 製品 \(p\) の \(1\) 単位を分解処理して生成される製品 \(q \in Parent_p\) の量（\(q\)-units）
\(RUB_r\): 資源 \(r\) の利用可能量上限（\(r\)-units）
\(R_{r}^{p}\): 製品 \(p\) の \(1\) 単位を（組み立て，分解，移動）処理する際に必要な資源 \(r\) の量（\(r\)-units）； ここで，\(r\)-unit とは，資源 \(r\) の \(1\) 単位を表す．
\(CT_{ijr}^p\): 枝 \((i,j)\) 上で資源 \(r\) を利用して製品 \(p\) を処理する際のサイクル時間（単位期間）
\(LT_{ijr}^p\): 枝 \((i,j)\) 上で資源 \(r\) を利用して製品 \(p\) を処理する際のリード時間（単位期間）
\(VAL_{i}^p\): 点 \(i\) 上での製品 \(p\) の価値（円）
\(SSR_i^p\): 点 \(i\) 上での製品 \(p\) の安全在庫係数．（無次元）
\(VAR_p\): 製品 \(p\) の変動比率（\(p\)-units）；これは，製品ごとの需要の分散と平均の比が一定と仮定したとき， 「需要の分散 \(/\) 需要の平均」と定義される値である．
\(ratio\): 利子率（% \(/\) 単位期間）
\(EIC_{ij}^p\): 枝 \((i,j)\) 上で資源 \(r\) を用いて処理（組み立て，分解，輸送）される 製品 \(p\) に対して定義されるエシェロン在庫費用（円 \(/\)単位期間）； この値は，以下のように計算される． \[ EIC_{ijr}^p =\max\{ ratio \times \left(VAL_{j}^p- \sum_{q \in Child_p} \phi_{qp} VAL_{i}^q \right)/ 100, 0 \} \ \ \ (i,j,r,p) \in Assemble \] \[ EIC_{ijr}^p =\max\{ ratio \times \left(\sum_{q \in Parent_p} \phi_{pq} VAL_{j}^q -VAL_{i}^p \right)/ 100, 0 \} \ \ \ (i,j,r,p) \in Disassemble \] \[ EIC_{ijr}^p =\max\{ ratio \times \left(VAL_{j}^p -VAL_{i}^p\right)/ 100, 0 \} \ \ \ (i,j,r,p) \in Transit \]
\(CFP_{ijr}\): 枝 \((i,j)\) で資源 \(r\) を使用したときの\(CO_2\)排出量 （g)； 輸送の場合には，資源 \(r\) の排出原単位 (g/km) に，枝 \((i,j)\) の距離 (km) を乗 じて計算しておく．
\(CFPV_{ijr}\): 資源 \(r\) の使用量（輸送の場合には積載重量）あたりの\(CO_2\)排出原単位（g \(/\) \(r\)-units)
\(CFPUB\): \(CO_2\)排出量上限（g）
変数は実数変数と\(0\)-\(1\)整数変数を用いる．
まず，実数変数を以下に示す．
\(w_{ijr}^p (\geq 0)\): 枝 \((i,j)\) で資源 \(r\) を利用して製品 \(p\) を処理する量を表す実数変数（\(p\)-units \(/\) 単位期間）
\(v_{ip}^+ (\geq 0)\): 点 \(i\) における製品 \(p\) の需要の超過量（需要が負のときには供給の超過量） を表す実数変数（\(p\)-units \(/\) 単位期間）
\(v_{ip}^- (\geq 0)\): 点 \(i\) における製品 \(p\) の需要の不足量（需要が負のときには供給の不足量） を表す実数変数（\(p\)-units \(/\) 単位期間）
\(0\)-\(1\)整数変数は，枝および枝上の資源の利用の有無を表現する．
定式化：
\[\begin{eqnarray*} \mbox{最小化} & \mbox{枝固定費用}+\mbox{枝・資源固定費用} + \\ & \mbox{枝・資源・製品変動費用}+ \mbox{供給量超過費用}+ \\ & \mbox{供給量不足費用}+ \mbox{需要量超過費用}+\mbox{需要量不足費用}+ \\ & \mbox{サイクル在庫費用} +\mbox{安全在庫費用}+\mbox{需要逸脱ペナルティ} \\ \mbox{条件} & \mbox{フロー整合条件} \\ & \mbox{資源使用量上限} \\ & \mbox{枝と枝上の資源の繋ぎ条件} \\ & \mbox{$CO_2$排出量上限制約} \end{eqnarray*}\]
\[ \mbox{枝・資源・製品変動費用}= \sum_{(i,j,r,p) \in ArcResProd} ARPVC_{ijr}^p w_{ijr}^p \]
\[ \mbox{需要量超過費用}= \sum_{(i,p) \in NodeProd %: D_{i}^{p}>0 } DPENALTY_{ip}^+ v_{ip}^+ \] \[ \mbox{需要量不足費用}= \sum_{(i,p) \in NodeProd %: D_{i}^{p}>0 } DPENALTY_{ip}^- v_{ip}^- \] \[ \mbox{サイクル在庫費用} = \sum_{(i,j,r,p) \in ArcResProd} \frac{EIC_{ijr}^p CT_{ijr}^p }{2} w_{ijr}^p \] \[ \mbox{安全在庫費用} = \sum_{(i,j,r,p) \in ArcResProd} ratio \times VAL_j^p SSR_i^p \sqrt{VAR_p LT_{ijr}^p w_{ijr}^p} \] 上式における平方根は区分的線形関数で近似するものとする． \[ \mbox{需要逸脱ペナルティ} = \sum_{(i,p) \in NodeProd} DPENALTY_{ip}^{+} v_{ip}^+ + DPENALTY_{ip}^{-} v_{ip}^- \]
一般化フロー整合条件 \[ \sum_{r \in Res, j \in N: (j,i,r,p) \in Transit \cup Assemble} w_{jir}^p+ \sum_{r \in Res, j \in N: (j,i,r,q) \in Disassemble, p \in Parent_q} \phi_{qp} w_{jir}^q \\ = \sum_{r \in Res, k \in N: (i,k,r,p) \in Transit \cup Disassemble} w_{ikr}^p+ \sum_{r \in Res, k \in N: (i,k,r,q) \in Assemble, p \in Child_q} \phi_{pq} w_{ikr}^q+ \\ (\mbox{ if } (i,p) \in NodeProd \mbox{ then } D_i^p -v_{ip}^{-} +v_{ip}^{+} \mbox{ else } 0) \ \ \ \forall i \in N, p \in Prod \]
資源使用量上限 \[ \sum_{p \in Prod: (i,j,r,p) \in ArcResProd} R_{r}^p w_{ijr}^p \leq RUB_{r} z_{ijr} \ \ \ \forall (i,j,r) \in ArcRes \]
枝と枝上の資源の繋ぎ条件 \[ z_{ijr} \leq y_{ij} \ \ \ \forall (i,j,r) \in ArcRes \]
\(CO_2\)排出量上限制約
\[ \sum_{(i,j,r) \in ArcRes} CFP_{ijr} z_{ijr} + \sum_{ (i,j,r,p) \in ArcResProd} CFPV_{ijr} R_r^p w_{ijr}^p \leq CFPUB \]
NodeData=[ "source1", "source2", "plantin", "plantout", "customer" ]
#Arc list, fixed cost (ArcFC) and Diatance (km)
ArcData={("source1", "plantin"): [100,200],
("source2", "plantin"): [100,500],
("plantin", "plantout"): [100,0],
("plantout", "customer"): [100,20],
}
#Prod data
#weight (ton), variability of product (VAR)= variance/mean
ProdData={
"apple": [0.01,0],
"melon": [0.012,0],
"bottle":[0.03,0],
"juice1":[0.05,2],
"juice2":[0.055,3],
}
#resource list, fixed cost（未使用）, upper bound
#corbon foot print (CFP) data (kg/km), variable term of corbon foot print (CFPV) (kg/ton km)
ResourceData={
"line1": [100,100,0, 0],
"line2": [20,100,0, 0],
"vehicle": [100,100,0.8, 0.1],
"ship": [30,180, 0.2, 0.1]
}
#Resource-Prod data (resource usage)
ResourceProdData = {
("line1", "juice1"): 1,
("line2", "juice1"): 1,
("line1", "juice2"): 1,
("line2", "juice2"): 1,
("vehicle", "juice1"): 1,
("ship", "juice1"): 1,
("vehicle", "juice2"): 1,
("ship", "juice2"): 1,
("vehicle", "apple"): 1,
("vehicle", "melon"): 1,
("ship", "apple"): 1,
("ship", "bottle"): 1,
}
#ArcResource list, fixed cost (ArcResourceFC)
ArcResourceData={
("source1", "plantin","vehicle"): 10,
("source2", "plantin","ship"): 30,
("plantin", "plantout","line1"): 50,
("plantin", "plantout","line2"): 100,
("plantout", "customer","vehicle"): 20,
("plantout", "customer","ship"): 40,
}
#ArcResourceProd data,
# type: 0=transport, 1=assemble, 2=dis-assemble,
# variable cost, cycle time, lead time (LT), and upper bound of flow volume (UB)
ArcResourceProdData= {
("source1", "plantin","vehicle","apple"): [0,1,1,10,50],
("source1", "plantin","vehicle","melon"): [0,2,2,20,50],
("source2", "plantin","ship","apple"): [0,3,3,15,50],
("source2", "plantin","ship","bottle"): [0,3,3,15,50],
("plantin", "plantout","line1","juice1"): [1,10,5,1,10],
("plantin", "plantout","line1","juice2"): [1,10,5,1,10],
("plantin", "plantout","line2","juice1"): [1,5,3,1,10],
("plantin", "plantout","line2","juice2"): [1,5,3,1,10],
("plantout", "customer","vehicle","juice1"): [0,10,5,3,10],
("plantout", "customer","vehicle","juice2"): [0,10,5,4,10],
("plantout", "customer","ship","juice1"): [0,2,2,10,10],
("plantout", "customer","ship","juice2"): [0,2,2,12,10],
}
#Node-Prod data; value, demand (negative values are supply), DPENALTY+, DPENALTY- (demand violation penalties)
NodeProdData ={
("source1", "apple"): [10,-30,0,10000],
("source1", "melon"): [10,-50,0,10000],
("source2", "apple"): [5,-100,0,10000],
("source2", "bottle"): [5,-100,0,10000],
("plantin", "apple"): [15,0,0,0],
("plantin", "melon"): [15,0,0,0],
("plantin", "bottle"): [8,0,0,0],
("plantout", "juice1"): [150,0,0,0],
("plantout", "juice2"): [160,0,0,0],
("customer", "juice1"): [170,10,0,10000],
("customer", "juice2"): [180,10,0,10000],
}
#BOM
Unit = {
("juice1", "apple"): 2,
("juice1", "melon"): 1,
("juice1", "bottle"): 1,
("juice2", "apple"): 2,
("juice2", "melon"): 2,
("juice2", "bottle"): 1
}
phi = {}
for (q,p) in Unit:
phi[p,q] = Unit[q,p]
ratio = 5.0 #interest ratio
SSR = 1.65 #safety stock ratio
CFPUB =400.0 #upper bound of carbon foot print (kg)
ArcList, ArcFC, Distance= gp.multidict(ArcData)
Prod, Weight, VAR = gp.multidict(ProdData)
Child, Parent ={}, {}
for (p,q) in phi:
if q in Child:
Child[q].append(p)
else:
Child[q] = [p]
if p in Parent:
Parent[p].append(q)
else:
Parent[p] = [q]
Res, ResourceFC, ResourceUB, CFP, CFPV = gp.multidict(ResourceData)
ResourceProd, R = gp.multidict(ResourceProdData)
ArcResource, ArcResourceFC = gp.multidict(ArcResourceData)
ArcResourcePair= gp.tuplelist([(i,j,r) for (i,j,r) in ArcResource])
ArcResourceProd, Type, VariableCost, CycleTime, LT, UB = gp.multidict(ArcResourceProdData)
ArcResourceProdPair= gp.tuplelist(ArcResourceProd)
TransPair, AsmblPair, DisasmblPair =[],[],[]
for (i,j,r,p) in Type:
if Type[i,j,r,p]==1:
AsmblPair.append( (i,j,r,p) )
elif Type[i,j,r,p]==2:
DisasmblPair.append( (i,j,r,p) )
else:
TransPair.append( (i,j,r,p) )
NodeProd, VAL, Demand, DP_plus, DP_minus =gp.multidict(NodeProdData)
NodeProdPair=gp.tuplelist(NodeProd)
DemandNodeProdPair=[(i,p) for (i,p) in Demand if Demand[i,p]>0]
SupplyNodeProdPair=[(i,p) for (i,p) in Demand if Demand[i,p]<0]products, arcs, nodes = {}, {}, {}
activities, modes, resources = {}, {}, {}
model = Model(name="Generic LND + Safety Stock Allocation")
for (i,j) in ArcList:
nodes[i] = Node(name=f"node({i})")
nodes[j] = Node(name=f"node({j})")
arcs[i,j] = model.addArc(name=f"arc({i},{j})", source= nodes[i], sink=nodes[j]) #need fixed_cost?
for (i,j,r) in ArcResource:
resources[i,j,r] = model.addResource(name=f"resource({i},{j},{r})", capacity=ResourceUB[r]) #need fixed_cost?
for p in Prod:
products[p] = model.addProduct(name=f"product({p})")
for (i,j,r,p) in ArcResourceProd:
if (i,j,p) not in activities:
if Type[i,j,r,p]==0: #輸送
activities[i,j,p] = model.addActivity(name=f"act({i},{j},{p})", atype="transport", product=products[p])
else: #製造
activities[i,j,p] = model.addActivity(name=f"act({i},{j},{p})", atype="make", product=products[p])
mode = Mode(name=f"mode({r})", variable_cost=VariableCost[i,j,r,p], fixed_cost= ArcResourceFC[i,j,r])
mode.addResource(resource = resources[i,j,r], variable= R[r,p])
if Type[i,j,r,p] == 1: #assemble
for q in Child[p]:
mode.addComponent(component = products[q], quantity= phi[q,p])
elif Type[i,j,r,p] == 2: #disassemble
for q in Parent[p]:
mode.addByproduct(byproduct = products[q], quantity= phi[p,q])
activities[i,j,p].addMode(mode)
arcs[i,j].addActivity(activity=activities[i,j,p])
for(i,p) in Demand:
if Demand[i,p]>0:
model.addData(dtype="demand", product=products[p], node=nodes[i], amount=Demand[i,p] )
elif Demand[i,p]<0: #供給
model.addData(dtype="supply", product=products[p], node=nodes[i], amount= -Demand[i,p], under_penalty =0. )#多期間モデル（在庫も考慮）に拡張
for i,t in enumerate(pd.date_range(start="2024/1/1", end="2024/2/28", freq="w")):
model.addPeriod(name= f"period({i})", start=t)
model.update()
planning_horizon = len(model.period_list)
print("T=",planning_horizon)
#工場の入口にボトルの在庫を置く
bottle = model.products["product(bottle)"]
plantin = model.nodes["node(plantin)"]
act_inv_bottle = model.addActivity(name="act bottle inventory", atype="inventory", product= bottle)
#mode_inv_bottle = Mode(name="mode bottle inventory") #在庫が資源を必要としたり，上限がある場合にはモードで定義
#act_inv_bottle.addMode(mode_inv_bottle)
plantin.addActivity(act_inv_bottle)
model.addData(dtype="value",product=bottle, node=plantin, amount=10.) #製品の価値を入力
model.addData(dtype="inventory",product=bottle, node=plantin, period= model.periods["period(0)"],amount=30.) #製品の在庫量を入力
model.addData(dtype="inventory",product=bottle, node=plantin, period= model.periods["period(7)"],amount=10.)
#計画期間を与えて， 過去の実績を定数として入力し，計画期間内だけ最適化
#工場の出口にジュース1の在庫を置く
juice1 = model.products["product(juice1)"]
plantout = model.nodes["node(plantout)"]
act_inv_juice1 = model.addActivity(name="act juice inventory", atype="inventory", product= juice1)
plantout.addActivity(act_inv_juice1)
model.addData(dtype="value",product=juice1, node=plantout, amount=40.)
model.interest_rate = 0.15 #在庫保管比率を設定
#model.period_index
#単一期間の場合にはstart=0, end=1を指定する．
start_period_name = 'period(1)'
end_period_name = 'period(7)'
gp_model = optimize(model, start_period_name, end_period_name)T= 8
計画期間 = 1 8
Welcome to the CBC MILP Solver
Version: 2.10.3
Build Date: Dec 15 2019
command line - /Users/mikiokubo/Documents/dev/scmopt2/env/lib/python3.10/site-packages/pulp/solverdir/cbc/osx/64/cbc /var/folders/c0/1l3dlhys02nbkv4fyj9pkm680000gn/T/adec1d2053e149c2ba4735a0c318a047-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/c0/1l3dlhys02nbkv4fyj9pkm680000gn/T/adec1d2053e149c2ba4735a0c318a047-pulp.sol (default strategy 1)
At line 2 NAME MODEL
At line 3 ROWS
At line 208 COLUMNS
At line 1325 RHS
At line 1529 BOUNDS
At line 1837 ENDATA
Problem MODEL has 203 rows, 307 columns and 627 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 2490.76 - 0.00 seconds
Cgl0003I 0 fixed, 0 tightened bounds, 42 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 1 strengthened rows, 0 substitutions
Cgl0004I processed model has 160 rows, 237 columns (84 integer (84 of which binary)) and 511 elements
Cbc0038I Initial state - 48 integers unsatisfied sum - 6.59365
Cbc0038I Pass 1: suminf. 0.34286 (2) obj. 1.50001e+08 iterations 56
Cbc0038I Solution found of 1.50001e+08
Cbc0038I Relaxing continuous gives 1.50001e+08
Cbc0038I Before mini branch and bound, 35 integers at bound fixed and 70 continuous
Cbc0038I Full problem 160 rows 237 columns, reduced to 80 rows 103 columns
Cbc0038I Mini branch and bound improved solution from 1.50001e+08 to 4555 (0.02 seconds)
Cbc0038I Round again with cutoff of 4360.55
Cbc0038I Pass 2: suminf. 5.66031 (40) obj. 4360.55 iterations 51
Cbc0038I Pass 3: suminf. 5.50475 (39) obj. 4360.55 iterations 1
Cbc0038I Pass 4: suminf. 2.40407 (21) obj. 4360.55 iterations 45
Cbc0038I Pass 5: suminf. 2.38733 (14) obj. 4360.55 iterations 37
Cbc0038I Pass 6: suminf. 1.55124 (10) obj. 4360.55 iterations 12
Cbc0038I Pass 7: suminf. 1.27950 (10) obj. 4360.55 iterations 1
Cbc0038I Pass 8: suminf. 1.15876 (8) obj. 4360.55 iterations 6
Cbc0038I Pass 9: suminf. 1.15876 (8) obj. 4360.55 iterations 0
Cbc0038I Pass 10: suminf. 1.15876 (8) obj. 4360.55 iterations 2
Cbc0038I Pass 11: suminf. 1.15876 (8) obj. 4360.55 iterations 0
Cbc0038I Pass 12: suminf. 1.15876 (8) obj. 4360.55 iterations 2
Cbc0038I Pass 13: suminf. 3.60267 (21) obj. 4360.55 iterations 40
Cbc0038I Pass 14: suminf. 3.30968 (21) obj. 4360.55 iterations 2
Cbc0038I Pass 15: suminf. 2.13019 (15) obj. 4360.55 iterations 16
Cbc0038I Pass 16: suminf. 2.13019 (15) obj. 4360.55 iterations 0
Cbc0038I Pass 17: suminf. 1.80162 (12) obj. 4360.55 iterations 6
Cbc0038I Pass 18: suminf. 1.80162 (12) obj. 4360.55 iterations 0
Cbc0038I Pass 19: suminf. 1.80162 (12) obj. 4360.55 iterations 2
Cbc0038I Pass 20: suminf. 1.80162 (12) obj. 4360.55 iterations 0
Cbc0038I Pass 21: suminf. 2.00162 (13) obj. 4360.55 iterations 14
Cbc0038I Pass 22: suminf. 2.05552 (12) obj. 4360.55 iterations 3
Cbc0038I Pass 23: suminf. 1.58203 (12) obj. 4360.55 iterations 1
Cbc0038I Pass 24: suminf. 2.05552 (12) obj. 4360.55 iterations 1
Cbc0038I Pass 25: suminf. 2.90590 (16) obj. 4360.55 iterations 17
Cbc0038I Pass 26: suminf. 2.90590 (16) obj. 4360.55 iterations 0
Cbc0038I Pass 27: suminf. 1.88733 (13) obj. 4360.55 iterations 7
Cbc0038I Pass 28: suminf. 1.88733 (13) obj. 4360.55 iterations 0
Cbc0038I Pass 29: suminf. 1.88733 (13) obj. 4360.55 iterations 2
Cbc0038I Pass 30: suminf. 1.88733 (13) obj. 4360.55 iterations 0
Cbc0038I Pass 31: suminf. 1.88733 (13) obj. 4360.55 iterations 2
Cbc0038I No solution found this major pass
Cbc0038I Before mini branch and bound, 7 integers at bound fixed and 71 continuous
Cbc0038I Full problem 160 rows 237 columns, reduced to 106 rows 121 columns - 6 fixed gives 95, 114 - still too large
Cbc0038I Full problem 160 rows 237 columns, reduced to 57 rows 73 columns
Cbc0038I Mini branch and bound did not improve solution (0.04 seconds)
Cbc0038I After 0.04 seconds - Feasibility pump exiting with objective of 4555 - took 0.04 seconds
Cbc0012I Integer solution of 4555 found by feasibility pump after 0 iterations and 0 nodes (0.04 seconds)
Cbc0038I Full problem 160 rows 237 columns, reduced to 113 rows 163 columns - 1 fixed gives 113, 162 - still too large
Cbc0038I Full problem 160 rows 237 columns, reduced to 113 rows 162 columns - too large
Cbc0031I 72 added rows had average density of 32.902778
Cbc0013I At root node, 72 cuts changed objective from 2610.5238 to 4403.656 in 100 passes
Cbc0014I Cut generator 0 (Probing) - 17 row cuts average 2.0 elements, 0 column cuts (49 active) in 0.069 seconds - new frequency is -100
Cbc0014I Cut generator 1 (Gomory) - 1905 row cuts average 122.2 elements, 0 column cuts (0 active) in 0.028 seconds - new frequency is 1
Cbc0014I Cut generator 2 (Knapsack) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.012 seconds - new frequency is -100
Cbc0014I Cut generator 3 (Clique) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.001 seconds - new frequency is -100
Cbc0014I Cut generator 4 (MixedIntegerRounding2) - 370 row cuts average 9.0 elements, 0 column cuts (0 active) in 0.029 seconds - new frequency is 1
Cbc0014I Cut generator 5 (FlowCover) - 11 row cuts average 3.7 elements, 0 column cuts (0 active) in 0.024 seconds - new frequency is -100
Cbc0014I Cut generator 6 (TwoMirCuts) - 223 row cuts average 33.9 elements, 0 column cuts (0 active) in 0.006 seconds - new frequency is 1
Cbc0010I After 0 nodes, 1 on tree, 4555 best solution, best possible 4403.656 (0.39 seconds)
Cbc0038I Full problem 160 rows 237 columns, reduced to 95 rows 142 columns - 1 fixed gives 94, 141 - still too large
Cbc0038I Full problem 160 rows 237 columns, reduced to 94 rows 141 columns - too large
Cbc0016I Integer solution of 4545 found by strong branching after 6131 iterations and 63 nodes (0.63 seconds)
Cbc0016I Integer solution of 4530 found by strong branching after 6132 iterations and 64 nodes (0.63 seconds)
Cbc0012I Integer solution of 4527.5 found by rounding after 6293 iterations and 68 nodes (0.63 seconds)
Cbc0004I Integer solution of 4525 found after 6507 iterations and 77 nodes (0.64 seconds)
Cbc0038I Full problem 160 rows 237 columns, reduced to 71 rows 114 columns
Cbc0004I Integer solution of 4505 found after 6986 iterations and 94 nodes (0.68 seconds)
Cbc0012I Integer solution of 4495 found by rounding after 7743 iterations and 115 nodes (0.74 seconds)
Cbc0016I Integer solution of 4490 found by strong branching after 7752 iterations and 118 nodes (0.75 seconds)
Cbc0001I Search completed - best objective 4489.999999999959, took 8573 iterations and 134 nodes (0.78 seconds)
Cbc0032I Strong branching done 1774 times (15375 iterations), fathomed 23 nodes and fixed 103 variables
Cbc0035I Maximum depth 13, 212 variables fixed on reduced cost
Cuts at root node changed objective from 2610.52 to 4403.66
Probing was tried 100 times and created 17 cuts of which 49 were active after adding rounds of cuts (0.069 seconds)
Gomory was tried 380 times and created 2302 cuts of which 0 were active after adding rounds of cuts (0.049 seconds)
Knapsack was tried 100 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.012 seconds)
Clique was tried 100 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.001 seconds)
MixedIntegerRounding2 was tried 380 times and created 803 cuts of which 0 were active after adding rounds of cuts (0.055 seconds)
FlowCover was tried 100 times and created 11 cuts of which 0 were active after adding rounds of cuts (0.024 seconds)
TwoMirCuts was tried 380 times and created 317 cuts of which 0 were active after adding rounds of cuts (0.018 seconds)
ZeroHalf was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Result - Optimal solution found
Objective value: 4490.00000000
Enumerated nodes: 134
Total iterations: 8573
Time (CPU seconds): 0.77
Time (Wallclock seconds): 0.78
Option for printingOptions changed from normal to all
Total time (CPU seconds): 0.77 (Wallclock seconds): 0.79
Opt. Val= 4505.0集合：
パラメータ：
変数：
定式化
\[ \begin{array}{l l l} \max & \tau & \\ s.t. & u_j \leq \sum_{i=(g,q)} \frac{1}{R_{ij}} y_{ij} & \forall j=(f,p), q \in CHILD_p \\ & \sum_{j} y_{ij} \leq u_i + I_i & \forall i \in N \\ & d_i \tau \leq u_i +I_i & \forall i \in N \\ & \sum_{p \in PRODUCT_f} u_{fp} \leq C_f \tau & \forall f \in V \\ & 0 \leq u_i \leq UB_i & \forall i \in N \\ & y_{ij} \geq 0 & \forall (i,j) \in A \end{array} \]
最初の制約（生産量と入庫輸送量の関係）： 工場 \(f\) における製品 \(p\) の生産量 \(u_{j} (j=(f,p))\) は、その部品 \(q\) の輸送量以下でなくてはいけない。 ここで、輸送量は出発地点 \(i\) における子製品 \(q\) の単位で表現されているので、親製品の量は変換比率 \(R_{ij}\) で割らなければならない。
2番目の制約（生産量と出庫輸送量の関係）： 点 \(i\) から出る輸送量 \(y_{ij}\) は、生産量 \(u_i\) とパイプライン在庫量 \(I_i\) の和以下でなくてはならない。
3番目の制約（需要満足条件）： 生産量は途絶期間内の需要量以上でなければならない。
4番目の制約（工場の生産容量制約）： 工場 \(f\) で生産される総量は、その容量 \(C_f\) 以下でなければならない。
各点（工場と製品の組）が途絶したシナリオ（点 \(i\) の生産量上限 \(UB_i\) を \(0\) に設定） において上の問題を解き、そのときの目的関数値が、点の途絶時に品切れがおきない最大の期間（余裕生存時間）となる。
BOM=nx.DiGraph() #BOM: bill of materials
# 工場グラフ
# 1 => 0
# 2 => 0
Capacity = {0: 300, 1: 500, 2: 200 }
Products = {0: ['P4','P5'], 1:['P1','P3'], 2: ['P2','P3']}
BOM.add_weighted_edges_from([ ('P1','P4', 1), ('P2','P5',2),('P3','P5',1) ])
model = Model(name = "Risk Exposure Index")
products, activities, modes, nodes, arcs, resources = {}, {}, {}, {}, {}, {}
#製品
for p in [f"P{i}" for i in range(1,6)]:
products[p] = model.addProduct(name=p)
#ノードとアーク
for i in range(3):
nodes[i] = model.addNode(name=i)
for (i,j) in [(1,0),(2,0)]:
arcs[i,j] = model.addArc(name=f"arc({i},{j})", source=nodes[i], sink=nodes[j])
#点の上に資源と活動を配置
for i, prods in Products.items():
resources[i] = model.addResource(name=f"res({i})", capacity=Capacity[i] )
for p in prods:
activities[i,p] = model.addActivity(name=f"act({i},{p})", atype="inventory", product = products[p])
modes[i,p] = Mode(name=f"mode({i},{p})", upper_bound = 100.)
for q in BOM.predecessors(p):
modes[i,p].addComponent(component=products[q], quantity=BOM[q][p]["weight"] )
modes[i,p].addResource(resource=resources[i], variable=1.)
activities[i,p].addMode(mode=modes[i,p])
nodes[i].addActivity(activities[i,p])
#パイプライン在庫，需要，製品の価値データ
for i, prods in Products.items():
for p in prods:
model.addData(dtype="inventory", product=products[p], node=nodes[i], amount= 300.)
if p in ["P4", "P5"]:
model.addData(dtype="value", product=products[p], node=nodes[i], amount= 4.)
model.addData(dtype="demand", product=products[p], node=nodes[i], amount= 100.)
else:
model.addData(dtype="value", product=products[p], node=nodes[i], amount= 1.)点集合 \(N\) は在庫地点を表し，枝 \((i,j) \in A\) が存在するとき，在庫地点 \(i\) は在庫地点 \(j\) に補充 を行うことを表す．複数の在庫地点から補充を受ける点においては，補充された各々の品目を用いて 別の品目を生産すると考える．
点 \(i\) における品目の生産時間は，\(T_i\) 日である． \(T_i\) には各段階における生産時間の他に，待ち時間および輸送時間も含めて考える．
このとき点 \(j\) は，複数の在庫地点から補充を受けるので，点 \(j\) が品目を発注してから，すべての品目が揃うまで生産を開始できない．
点上で生産される品目は，点によって唯一に定まるものと仮定する． このとき，点と品目を表す添え字は同一であると考えられるので， 有向グラフ \(G=(N,A)\) は部品展開表を意味することになる． \((i,j) \in A\) のとき，点 \(j\) 上の品目 \(j\) は，点 \(i\) から補充される品目 \(i\) をもとに生産される． 品目 \(j\) を生産するのに必要な品目 \(i\) の量を \(\phi_{ij}\) と記す．
需要は，後続する点をもたない点 \(j\) 上で発生するものとし， その \(1\) 日あたりの需要量は，期待値 \(\mu_j\) の定常な分布をもつものとする． 点 \(j\) における \(t\) 日間における需要の最大値を \(D_j(t)\) 記す．
直接の需要をもたない点（後続点をもつ点） \(i\) に対する需要量の期待値 \(\mu_i\) は， \[ \mu_i = \sum_{(i,j) \in A} \phi_{ij} \mu_j \] と計算できる．
点 \(i\) に対する\(t\) 日間における需要の最大値 \(D_i(t)\) は， \[ D_i(t)= t \mu_i +\left( \sum_{(i,j) \in A} \phi_{ij} (D_j(t)-t \mu_j)^p \right)^{1/p} \] と計算されるものと仮定する．
ここで，\(p (\geq 1)\) は定数である． \(p\) が大きいほど，需要の逆相関が強いことを表す． 点 \(j\) 上での需要が常に同期していると仮定した場合には \(p=1\) が， 点 \(j\) 上での需要が独立な正規分布にしたがうと仮定した場合には \(p=2\) が推奨される．
点 \(i\) の保証リード時間を \(L_i\) と記す． ここでは，各点 \(i\) に対して，保証リード時間の下限は \(0\) とし，上限は \(\bar{L}_i\) であるとする．
点 \(j\) が品目を発注してから，すべての品目が 揃うまでの時間（日数）を入庫リード時間とよび， \(LI_j\) と記す．
点 \(j\) における入庫リード時間 \(LI_j\) は，以下の式を満たす． \[ L_i \leq LI_j \ \ \ \forall (i,j) \in A \]
入庫リード時間 \(LI_i\) に生産時間 \(T_i\) を加えたものが，補充の指示を行ってから生産を完了するまでの時間となる． これを，補充リード時間とよぶ． 補充リード時間から、保証リード時間 \(L_i\) を減じた時間内の最大需要に相当する在庫を保持していれば，在庫切れの心配がないことになる． 補充リード時間から \(L_i\) を減じた時間（\(L_{i+1}+T_i-L_i\)）を正味補充時間とよぶ．
点 \(i\) における安全在庫量 \(I_i\) は，正味補充時間内における最大需要量から平均需要量を減じた量であるので， \[ I_i= D(LI_{i}+T_i-L_i) - (LI_{i}+T_i-L_i) \mu_i \] となる．
記述を簡単にするために，点 \(i\) において，保証リード時間が \(L\)， 入庫リード時間が \(LI\) のときの安全在庫費用を表す，以下の関数 \(HC_i(L,LI)\) を導入しておく．
\[ HC_i (L,LI)=h_i \left\{ D(LI+T_i-L) - (LI+T_i-L) \mu_i \right\} \]
上の記号を用いると，木ネットワークモデルにおける安全在庫配置問題は，以下のように定式化できる．
\[ \begin{array}{l l l } minimize & \sum_{i \in N} HC_i (L_i,LI_i) & \\ s.t. & L_i \leq LI_{i}+T_i & \forall i \in N \\ & L_i \leq LI_j & \forall (i,j) \in A \\ & 0 \leq L_i \leq \bar{L}_i & \forall i \in N \end{array} \]
ここで，最初の制約式は，正味補充時間が非負であることを表す．
tabu_search_for_SSA (G, ProcTime, z, mu, sigma, h, LTUB, max_iter=100, TLLB=1, TLUB=10, seed=1)
一般のネットワークの安全在庫配置モデルに対するタブー探索（mu未使用；NRT日の最大在庫量をシミュレーションもしくは畳み込みによって事前計算？もしくは正規分布で近似）
#
n = 7
G = SCMGraph()
for i in range(n):
G.add_node(i)
G.add_edges_from([(0, 2), (1, 2), (2,4), (3,4), (4,5), (4,6)])
#点のラベルの付け替え
mapping ={i:idx for idx,i in enumerate(G)}
G = nx.relabel_nodes(G, mapping=mapping, copy=True)
z = np.full(len(G),1.65)
h = np.array([1,1,3,1,5,6,6])
mu = np.array([200,200,200,200,200,100,100])
sigma = np.array([14.1,14.1,14.1,14.1,14.1,10,10])
ProcTime = np.array([6,2,3,3,3,3,3], int)
LTUB = np.array([0,0,0,0,0,4,3],int) #最大保証リード時間（需要地点のみ意味がある）
max_iter = 10
seed =123
TLLB, TLUB = 2,10 # tabu length is a random number between (TLLB, TLUB)
best_cost, best_sol, best_NRT, best_MaxLI, best_MinLT = tabu_search_for_SSA(G, ProcTime, z, mu, sigma,
h, LTUB, max_iter = 10, TLLB =1, TLUB =3, seed = seed)
print("最良値", best_cost)
print("最良解", best_sol)
print("正味補充時間", best_NRT)
print("最大補充リード時間", best_MaxLI)
print("最小保証リード時間", best_MinLT)最良値 374.82595172371384
最良解 [1 1 0 0 1 0 0]
正味補充時間 [6. 2. 0. 0. 6. 0. 0.]
最大補充リード時間 [6. 2. 3. 3. 6. 3. 3.]
最小保証リード時間 [0. 0. 3. 3. 0. 4. 3.]起動停止問題(unit commitment problem)とは，発電機の起動と停止を最適化するためのモデルであり， 各日の電力需要を満たすように時間ごとの発電量を決定する．
定式化のためには，動的ロットサイズ決定問題に以下の条件を付加する必要がある．
この制約の定式化が実用化のための鍵になるので，簡単に説明しておく．
集合:
変数:
\(y_t^p\): 期 \(t\) に発電機 \(p\) が稼働するとき \(1\)
\(z_t^p\): 期 \(t\) に発電機 \(p\) が稼働を開始(switch-on)するとき \(1\)
\(w_t^p\): 期 \(t\) に発電機 \(p\) が停止を開始(switch-off)するとき \(1\)
上の変数の関係は以下の式で表すことができる．
\[ \begin{array}{l l } z_t^p \leq y_t^p & \forall p \in P, t=1,2, \ldots,T \\ z_t^p -w_{t}^p = y_t^p -y_{t-1}^p & \forall p \in P, t=1,2, \ldots,T \end{array} \]
\(\Rightarrow\) 「期 \(t\) に稼働していない」 ならば 「\(t-\alpha+1\) から \(t\) までは開始できない」
\(\Rightarrow\) \(y_t^p\) ならば \(z_{s}^p=0\) (\(\forall s=t-\alpha+1,\ldots,t\))
\[ \displaystyle\sum_{s=t-\alpha+1}^t z_{s}^p \leq y_t^p \ \ \ \forall t=\alpha,\alpha+1, \ldots,T \]
弱い定式化：
\(\Rightarrow\) 「期 \(t\) に開始した」 ならば 「\(t\) から\(t+\alpha+1\) までは稼働」
\[ \alpha z_{t}^p \leq \displaystyle\sum_{s=t}^{t+\alpha-1} y_t^p \ \ \ \forall t=1,2, \ldots,T+1-\alpha \]
\(\Rightarrow\) 「期 \(t\) に稼働した」 ならば 「\(t-\beta+1\) から \(t\) までは終了できない」
\(\Rightarrow\) \(y_t^p=1\) ならば \(w_{s}^p=0\) (\(\forall s=t-\beta+1,\ldots,t\))
\[ \displaystyle\sum_{s=t-\beta+1}^{t} w_{s}^p \leq 1-y_t^p \ \ \ \forall t=\beta,\beta+1, \ldots,T \]
弱い定式化：
\(\Rightarrow\) 「期 \(t\) に停止した」 ならば 「\(t\)から \(t+\beta-1\) までは稼働しない」
\[ \beta w_{t}^p \leq \displaystyle\sum_{s=t}^{t+\beta-1} (1-y_t^p) \ \ \ \forall t=1,2, \ldots,T+1-\beta \]