<jupyter_start><jupyter_text>Supply Chain Optimization Modelshttp://www.logopt.com/mikiokubo/lecture/lecture4.htm 準備<jupyter_code>from amplpy import AMPL, Environment, DataFrame,  register_magics
ampl = AMPL(Environment(r'/Applications/ampl_macosx64/'))
register_magics(store_name='_ampl_cells', ampl_object=ampl)<jupyter_output><empty_output><jupyter_text>ナップサック問題AMPLのコマンドやモデルを入力するには，セルの最初に```python%%ampl_evalreset;```を入れる．<jupyter_code>%%ampl_eval
# knapsack problem
reset;
reset data;
set Bear;
var x{Bear} binary; 
param Value{Bear};
param Weight{Bear};
param MaxWeight;
maximize Profit: sum{i in Bear} Value[i]*x[i];
subject to Capacity: sum{i in Bear} Weight[i]*x[i]<=MaxWeight;<jupyter_output><empty_output><jupyter_text>データは amplオブジェクトのsetやparam属性で代入する．<jupyter_code># knapsack problem data
# set Bear := VerySmall Small Medium Large;
# param: Value Weight := 
# VerySmall 16 2 
# Small     19 3
# Medium    23 4 
# Large     28 5; 
# param MaxWeight := 7; 

ampl.set['Bear'] = ["VerySmall","Small","Medium","Large"]
ampl.param['Value'] = [16, 19, 23, 28]
ampl.param['Weight'] = [2,3,4,5]
ampl.param['MaxWeight'] = 7
%%ampl_eval
display Bear;
display Value;
display Weight;
display MaxWeight;
%%ampl_eval
option solver gurobi;
solve;
%%ampl_eval
show;
display Profit;
display x;
expand Profit;
expand Capacity;
ampl.exportData("kp.dat")<jupyter_output><empty_output><jupyter_text>巡回セールスマン問題<jupyter_code>ampl = AMPL(Environment(r'/Applications/ampl_macosx64/'))
register_magics(store_name='_ampl_cells', ampl_object=ampl)
%%ampl_eval
# MTZ formulation for asymmetric traveling salesman problem
reset;
param n;
set V :=  1..n;   
set V0 := 2..n;
set A :=V cross V;  
param D { A } >= 0, default 9999999;  
var x { A } binary ;   
var u { V0 } >=1,<=n-1;
minimize total_cost:
   sum {(i,j) in A} D[i,j] * x[i,j]; 
subject to Degree1 {i in V}:
   sum {(i,j) in A } x[i,j] =1 ;
subject to Degree2 {i in V}:
   sum {(j,i) in A } x[j,i] =1 ;
subject to MTZ{ (i,j) in A: i != j and j!=1 and i!=1}:
   u[i]+1 -(n-1)*(1-x[i,j])<=u[j]; 
%%ampl_eval
# MTZ formulation for asymmetric traveling salesman problem
# lifted MTZ constraints
reset;
param n >=0;
set V :=  1..n ;   
set V0 := 2..n;
set A :=V cross V;  
param D { A } >= 0, default 9999999;  
var x { A } binary ;   
var u { V0 } >=1,<=n-1;
minimize total_cost:
   sum {(i,j) in A} D[i,j] * x[i,j]; 
subject to Degree1 {i in V}:
   sum {(i,j) in A } x[i,j] =1 ;
subject to Degree2 {i in V}:
   sum {(j,i) in A } x[j,i] =1 ;
subject to MTZ{ (i,j) in A: i != j and j!=1 and i!=1}:
   u[i]+1 -(n-1)*(1-x[i,j]) + (n-3)*x[j,i]<=u[j]; 
subject to LiftedLB{ i in V0}:
   1+(1-x[1,i]) +(n-3)*x[i,1] <= u[i];
subject to LiftedUB{ i in V0}:
 u[i] <=(n-1)-(1-x[i,1])-(n-3)*x[1,i]; 
import random
n = 5
ampl.param['n'] = n
D = ampl.getParameter("D")
for i in range(1,n+1):
    for j in range(1,n+1):
        if i!=j:
            D[i,j] = random.randint(1,10)
%%ampl_eval
display n;
display D;
%%ampl_eval
expand total_cost;
option solver gurobi;
solve;
%%ampl_eval
display x;<jupyter_output>x :=
1 2   1
2 4   1
3 5   1
4 3   1
5 1   1
;<jupyter_text>施設配置問題<jupyter_code>ampl = AMPL(Environment(r'/Applications/ampl_macosx64/'))
register_magics(store_name='_ampl_cells', ampl_object=ampl)
%%ampl_eval
reset;
set Warehouse;   
set Market;  
set Edge := Warehouse cross Market;
param Supply {Warehouse} >= 0, default 9999999;  
param Demand {Market} >= 0;  
param TransCost { Edge } >= 0; 
param FixedCost {Warehouse};
var x { Edge }  >= 0;   
var y {Warehouse } binary;
minimize total_cost:
   sum {(i,j) in  Edge} TransCost[i,j] * x[i,j] + sum{ i in Warehouse } FixedCost[i]* y[i];
subject to SupplyUpperBound {i in Warehouse}:
   sum {j in Market} x[i,j] <= Supply[i]*y[i];
subject to SupplyUpperBound2 {(i,j) in Edge}:
   x[i,j] <= Demand[j]*y[i];
subject to DemandConst {j in Market}:
   sum {i in Warehouse} x[i,j] = Demand[j];
# set Warehouse :=A B C;
# set Market :=  N  K  O   S  T; 
# param Demand := N 180 K 80 O 200 S 160 T 220;
# param TransCost: N K O S T := 
#      A 1000  800   600  500  400   
#      B 600   500   400  300  600   
#      C 300   400   500  500  900 ;
# param FixedCost := A 10000 B 10000 C 10000;  
tc = [ [1000,  800,   600,  500,  400],   
     [ 600,   500,   400,  300,  600],   
     [ 300,   400,   500,  500,  900] ]
warehouse = ["A","B","C"]
market = ["N","K","O","S","T"]
ampl.set['Warehouse'] = warehouse
ampl.set['Market'] = market 
ampl.param['Demand'] = [180, 80, 200, 160, 220]
ampl.param['FixedCost'] = [10000, 10000, 10000]
ampl.param['TransCost'] = {(w,m): tc[i][j] for i,w in enumerate(warehouse) for j,m in enumerate(market)}
%%ampl_eval
display Demand;
display TransCost;
%%ampl_eval
option solver gurobi;
solve;
display x;
display y;<jupyter_output>Gurobi 9.0.2: optimal solution; objective 332000
2 simplex iterations
x :=
A T   220
B O   200
B S   160
C K    80
C N   180
;

y [*] :=
A  1
B  1
C  1
;<jupyter_text>ロットサイズ決定問題通常の定式化とエシェロン在庫モデル<jupyter_code>ampl = AMPL(Environment(r'/Applications/ampl_macosx64/'))
register_magics(store_name='_ampl_cells', ampl_object=ampl)
%%ampl_eval
# Big Bucket Lot Sizing Problem
reset;

param T;                     #number of periods 
set Period :={1..T}; 
set Prod;
set Parent {Prod} default {};
set Resource;
set ResourceProdPair within {Resource, Prod};

param SetUpCost{ Prod, Period };
param SetUpTime{ Prod, Period};
param VariableCost{ Prod, Period };
param HoldingCost{ Prod, Period };
param Demand{Prod,Period} default 0;
set DemandProd := setof {p in Prod, t in Period: Demand[p,t]>0} p;
param MinResourceUB{Prod, Period }>=0;
param Unit{p in Prod, Parent[p]};
param ResourceUB{ Resource, Period }>=0;
param R {ResourceProdPair} >=0;
param DemandSlackPenalty{DemandProd} default 99999;
param DemandSurplusPenalty{DemandProd} default 99999;
param FirstInv{Prod} default 0;
param LastInv{Prod} default 0;

var y{Prod,Period} binary; 
var x{Prod,Period} >=0;
var inv{Prod,Period diff {T} } >=0;
var vminus{DemandProd,Period}>=0;
var vplus{DemandProd,Period}>=0;

var TotalCost{2..5};

subject to ProductionSetupConnect{ p in Prod, t in Period }:
  x[p,t] <= MinResourceUB[p,t]*y[p,t]  ;

minimize cost:
sum{p in Prod,t in Period} SetUpCost[p,t]*y[p,t]+   sum{i in {2..5}} TotalCost[i];

#subject to Cost1:
#  TotalCost[1] = sum{p in Prod,t in Period} SetUpCost[p,t]*y[p,t] ;

subject to Cost2:
  TotalCost[2] = sum{p in Prod,t in Period}  VariableCost[p,t]*x[p,t];

subject to Cost3:
  TotalCost[3] = sum{p in Prod,t in Period diff {T} } HoldingCost[p,t]*inv[p,t];

subject to Cost4:
  TotalCost[4] = sum{p in Prod,t in Period}  DemandSlackPenalty[p]*vminus[p,t];

subject to Cost5:
  TotalCost[5] = sum{p in Prod,t in Period} DemandSurplusPenalty[p]*vplus[p,t];

subject to FlowConservation{p in Prod,t in Period}:
(if t=1 then FirstInv[p] else inv[p,t-1])+
    x[p,t] =
(if t=T then LastInv[p] else inv[p,t]) +
+sum{q in Parent[p]} Unit[p,q]*x[q,t]
+ (if p in DemandProd then 
   Demand[p,t]+vminus[p,t]-vplus[p,t] else 0) 
;

subject to ResourceUpperBound {r in Resource, t in Period}:
sum {p in Prod: (r,p) in ResourceProdPair}
   (R[r,p] * x[p,t] + SetUpTime[p,t]*y[p,t]) <= ResourceUB[r,t];

subject to Cut{p in Prod, t in Period diff {T} }:
  x[p,t] <=Demand[p,t]*y[p,t]+inv[p,t];

ampl.readData("lot1.dat")
%%ampl_eval
option solver gurobi;
solve;
display x;
ampl = AMPL(Environment(r'/Applications/ampl_macosx64/'))
register_magics(store_name='_ampl_cells', ampl_object=ampl)
%%ampl_eval
reset;
# Big Bucket Lot Sizing Problem
# Echelon Inventory Formulation
param T;                     #number of periods 
set Period :={1..T}; 
set Prod;
set Parent {Prod} default {};
set Resource;
set ResourceProdPair within {Resource, Prod};
param SetUpCost{ Prod, Period };
param SetUpTime{ Prod, Period};
param VariableCost{ Prod, Period };
param HoldingCost{ Prod, Period }; #Echelon Holding Cost
#Echelon Demand (sum of direct demand plus demand of ancestors)
param Demand{Prod,Period} default 0; 
set DemandProd := setof {p in Prod, t in Period: Demand[p,t]>0} p;

param MinResourceUB{Prod, Period }>=0;

param Unit{p in Prod, Parent[p]};

param ResourceUB{ Resource, Period }>=0;


param R {ResourceProdPair} >=0;
param DemandSlackPenalty{DemandProd} default 9999999;
param DemandSurplusPenalty{DemandProd} default 9999999;

param FirstInv{Prod} default 0;
param LastInv{Prod} default 0;
var y{Prod,Period} binary; 
var x{Prod,Period} >=0;
var inv{Prod,Period} >=0; #Echelon Inventory
var vminus{DemandProd,Period}>=0;
var vplus{DemandProd,Period}>=0;

#var TotalSetUpCost >=0;
#var TotalVariableCost >=0;
#var TotalHoldingCost>=0;
#var TotalSlack>=0;
#var TotalSurplus>=0;
var TotalCost{1..5};


#subject to Cost1:
#  TotalCost[1] = sum{p in Prod,t in Period} SetUpCost[p,t]*y[p,t] ;

subject to Cost2:
  TotalCost[2] = sum{p in Prod,t in Period}  VariableCost[p,t]*x[p,t];

subject to Cost3:
  TotalCost[3] = sum{p in Prod,t in Period} HoldingCost[p,t]*inv[p,t];

subject to Cost4:
  TotalCost[4] = sum{p in Prod,t in Period}  DemandSlackPenalty[p]*vminus[p,t];

subject to Cost5:
  TotalCost[5] = sum{p in Prod,t in Period} DemandSurplusPenalty[p]*vplus[p,t];


minimize cost:
#   TotalSetUpCost+TotalVariableCost+TotalHoldingCost+TotalSlack+TotalSurplus;
 sum{p in Prod,t in Period} SetUpCost[p,t]*y[p,t] +  sum{i in {2..5}} TotalCost[i];

subject to ProductionSetupConnect{ p in Prod, t in Period }:
  x[p,t] <= MinResourceUB[p,t]*y[p,t]  ;

subject to FlowConservation{p in Prod,t in Period}:
(if t=1 then FirstInv[p] else inv[p,t-1])+
    x[p,t] 
 =
(if t=T then LastInv[p] else inv[p,t]) +
 (if p in DemandProd then 
   Demand[p,t]+vminus[p,t]-vplus[p,t] else 0) 
;

#Echelon Connection Constraints
subject to EchelonConnect{p in Prod, t in Period}:
inv[p,t] >=sum{q in Parent[p]} Unit[p,q]*inv[q,t];

subject to ResourceUpperBound {r in Resource, t in Period}:
sum {p in Prod: (r,p) in ResourceProdPair}
   (R[r,p] * x[p,t] + SetUpTime[p,t]*y[p,t]) <= ResourceUB[r,t];

subject to Cut{p in Prod, t in Period}:
  x[p,t] <=Demand[p,t]*y[p,t]+inv[p,t];

ampl.readData("lot1.dat")
%%ampl_eval
option solver gurobi;
solve;
display x;
%%ampl_eval
### MELOS-R (Reverse)  AMPL File ### 
reset;
### Sets ### 
set Prod; 
set Parent{Prod} within Prod;
set Child{Prod}  within Prod;
set Cust 'Customers'; #Set of customers
set DC 'Warehouses';  #Set of warehouses
set Plnt 'Plants';    #Set of plants
set AllLine;
set Sply 'Suppliers'; #Set of suppliers
set RevSply;
set Node := Cust union DC union Plnt union Sply union RevSply;
set Route within {Node,Node,Prod};

set Arc  := setof { (i,j,p) in Route} (i,j); 
set Line{Plnt} within AllLine;
set ProdInLine{AllLine} within Prod;
set Pattern := 
setof{i in Plnt, ell in Line[i], p in ProdInLine[ell] } (i,ell,p);

set ProdIn{i in Plnt } :=
    setof { (j,i,p) in Route} p; 
set ProdOut{i in Plnt} :=
    setof {(i,j,p) in Route} p;

### Parameters ###
# Product data
param VOL{ Prod } >=0;
param HEIGHT{ Prod } >=0;
param Unit{p in Prod, Child[p] } >=0;
param RevUnit{ p in Prod, Parent[p]} >=0;

# Supplier data
param SupplyLB {i in Sply,p in Prod} >=0;
param SupplyUB {i in Sply,p in Prod} >=SupplyLB[i,p];
param SC {Sply,Prod} >=0;

# Reverse Supplier data
param RevSupplyLB {i in RevSply,p in Prod} >=0;
param RevSupplyUB {i in RevSply,p in Prod} >=RevSupplyLB[i,p];
param RevSC {RevSply,Prod} >=0;

# Warehouse data
 param FixedCost {DC union Plnt union AllLine} >=0;
 param yLB{DC union Plnt}>=0,<=1;
 param yUB{j in DC union Plnt}>=yLB[j], <=1;
 param VariableCost {DC,Prod} >=0;
 param DCLB {DC} >=0;
 param DCUB {j in DC} >= DCLB[j];
 param VAL{DC,Prod} >=0;
 param SAFETY{DC,Prod} >=0;

# Line data
param LLB {AllLine};
param LUB {ell in AllLine} >=LLB[ell];
param PC{ell in AllLine, ProdInLine[ell] };
param PLB{ell in AllLine, ProdInLine[ell] };
param PUB{ell in AllLine, p in ProdInLine[ell] }>= PLB[ell,p];
param R{ell in AllLine, ProdInLine[ell] };
param ITRL{ell in AllLine, ProdInLine[ell] } >0;
param LHoldingCost{ell in AllLine, ProdInLine[ell] }>=0;
param zLB{i in Plnt,ell in Line[i]} >=0, <=1;
param zUB{i in Plnt,ell in Line[i]} >=zLB[i,ell], <=1;
param VALL {ell in AllLine, ProdInLine[ell] } >=0;
param SAFETYL {ell in AllLine, ProdInLine[ell]} >=0;

# Customer data
param Demand {Cust,Prod}>=0;
param BOC    {Cust,Prod}>=0;
param RevDemand {Cust,Prod}>=0;
param RevBOC {Cust,Prod} >=0;

# Transportation data
param TransCost {Route} >=0;
param ITR {Route} >0;
param HoldingCost {Route} >=0;

# Other parameters
param ratio;
param SSR;

# Induced Set
set FinalProd := setof{k in Cust,p in Prod: Demand[k,p] >0 } p;
set RevProd   := setof{k in Cust,p in Prod: RevDemand[k,p] >0 } p;
set InitialProd := setof{s in Sply,p in Prod: SupplyUB[s,p]>0} p;
set RevInitialProd := setof{s in RevSply,p in Prod: RevSupplyUB[s,p]>0} p;

### Variables ###
var y 'open or not' {j in DC union Plnt} binary, >=yLB[j], <=yUB[j] ; 
var z {i in Plnt, ell in Line[i]} binary, >=zLB[i,ell],<=zUB[i,ell];
var w {i in Plnt, ell in Line[i], p in ProdInLine[ell]} >=0;
var x  'flow volume' {Route} >=0;
var xi 'backorder'   {Cust,FinalProd} >=0;
var revxi 'rev backorder' {Cust,RevProd} >=0;

### Objective ###
minimize cost: 
     sum {i in Sply, p in InitialProd, (i,j,p) in Route} SC[i,p] * x[i,j,p] +
     sum {i in RevSply, p in RevInitialProd, (j,i,p) in Route} RevSC[i,p] * x[j,i,p] +
     sum {(i,j,p) in Route} TransCost[i,j,p] * x[i,j,p] +
     sum {j in (DC union Plnt)} FixedCost[j] * y[j] +
     sum {i in Plnt, ell in Line[i]} FixedCost[ell] * z[i,ell] +
     sum {j in DC,  (j,k,p) in Route} VariableCost[j,p] * x[j,k,p] +
     sum {i in Plnt, ell in Line[i], p in ProdInLine[ell]} PC[ell,p] *w[i,ell,p]+
     sum {i in Plnt, ell in Line[i], p in ProdInLine[ell]} LHoldingCost[ell,p]/ITRL[ell,p]* w[i,ell,p]+
     sum {i in Plnt, ell in Line[i], p in ProdInLine[ell]} ratio*VALL[ell,p]*SAFETYL[ell,p] * w[i,ell,p]+
     sum {(i,j,p) in Route} HoldingCost[i,j,p] /ITR[i,j,p] * x[i,j,p] + 
     sum {j in DC,  (j,k,p) in Route} ratio*VAL[j,p]*SAFETY[j,p] * x[j,k,p]+
     sum {k in Cust, p in FinalProd} BOC[k,p]*xi[k,p]+
     sum {k in Cust, p in RevProd} RevBOC[k,p]*revxi[k,p]  ;

### Constraints ###
DemSat 'all demands are satisfied' {k in Cust, p in FinalProd}:
     sum {(j,k,p) in Route} x[j,k,p] + xi[k,p] = Demand[k,p];

RevDemSat 'all reverse demands are satisfied' {k in Cust, p in RevProd}:
     sum {(k,j,p) in Route} x[k,j,p] + revxi[k,p] = RevDemand[k,p];

Supplybnds 'upper and lower bounds of supply' {i in Sply, p in InitialProd}:
    SupplyLB[i,p] <= sum { (i,j,p) in Route} x[i,j,p] <= SupplyUB[i,p];

RevSupplybnds 'upper and lower bounds of reverse supply' {i in RevSply, p in RevInitialProd}:
    RevSupplyLB[i,p] <= sum { (j,i,p) in Route} x[j,i,p] <= RevSupplyUB[i,p];

DCLowerBound 'lower bound of warehouse' {j in DC}:
    DCLB[j] * y[j] <= sum {k in Node,p in Prod: (j,k,p) in Route} 
    VOL[p] * SSR /HEIGHT[p] *(1/ITR[j,k,p] +SAFETY[j,p])* x[j,k,p];

DCUpperBound 'upper bound of warehouse' {j in DC}:
  sum {k in Node, p in Prod: (j,k,p) in Route}
    VOL[p] * SSR /HEIGHT[p] * x[j,k,p] <=DCUB[j] * y[j];

flowcons 'flow conservation equation' {j in DC, p in Prod}:
   sum { (i,j,p) in Route} x[i,j,p] = 
    sum {(j,k,p) in Route} x[j,k,p];

LineResourceLB {i in Plnt, ell in Line[i]}:
LLB[ell] * z[i,ell] <= 
sum{ p in ProdInLine[ell] } R[ell,p]*w[i,ell,p];

LineResourceUB {i in Plnt, ell in Line[i]}:
sum{ p in ProdInLine[ell] } R[ell,p]*w[i,ell,p] 
<= LUB[ell] * z[i,ell];

LineLB {i in Plnt, ell in Line[i], p in ProdInLine[ell]}:
PLB[ell,p]* z[i,ell] <= w[i,ell,p];

LineUB {i in Plnt, ell in Line[i], p in ProdInLine[ell]}:
 w[i,ell,p] <= PUB[ell,p]* z[i,ell];

ConnectLineAndPlnt {i in Plnt, ell in Line[i]}:
z[i,ell] <= y[i];

flowconsPlntIn {i in Plnt, q in ProdIn[i]}:
sum{(s,i,q) in Route} x[s,i,q] =
 sum{ell in Line[i]: (i,ell,q) in Pattern} w[i,ell,q] +
 sum{ell in Line[i], p in ProdInLine[ell]: q in Child[p]} Unit[p,q]*w[i,ell,p];

flowconsPlntout {i in Plnt, r in ProdOut[i]}:
sum{ell in Line[i]: (i,ell,r) in Pattern} w[i,ell,r] +
sum{ell in Line[i], p in ProdInLine[ell]: r in Parent[p]}
 RevUnit[p,r]*w[i,ell,p]
= sum{(i,j,r) in Route} x[i,j,r];<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup