MPライブラリ

MPライブラリの使用法

ampl_notebook

Google Colab.で使う場合

ampl_notebook関数で、使うソルバー(modules)とlicense_uuidを設定する。

# !pip install amplpy
# from amplpy import ampl_notebook
# ampl = ampl_notebook(
#     modules=["highs","gurobi","cbc","scip","coin","gecode"],  #coin includes ipopt, couenne, bonmin
#     license_uuid=None)

ローカルで amplpyを動かす方法

AMPLをインストールしてampl.exeがある場所をEnvironmentで指定する。

from amplpy import AMPL, Environment, tools
ampl = AMPL(Environment("/Users/mikiokubo/Documents/ampl/"))

マジックコマンド ampl_eval

マジックコマンド%%ampl_evalを入れると、そのままAMPLのモデルとコマンドが記述できる。

reset;
var x1 integer, >= 0;
var x2 integer,  >= 0;
maximize Z: 3*x1 + 2*x2;
subject to Constraint1: 2*x1 + x2 <= 10;
subject to Constraint2: x1==1 or x2 ==5;
option solver gecode;
solve;

display x1,x2,Z;
gecode 6.2.0: optimal solution
13 nodes, 1 fails, objective 19
x1 = 1
x2 = 8
Z = 19

記述したamplモデルは、Pythonのamplインスタンス（グローバル変数）ampl に保管される。

たとえば、evalメソッドで、amplのshow;コマンドを渡すと、モデルの概要が表示される。

ampl.eval("show;")

variables:   x1   x2

constraint:   Constraint1

logical constraint:   Constraint2

objective:   Z

MPサポート

拡張されたMPソルバーインターフェースライブラリは、以下の演算子および式カテゴリーに対する新たなサポートを提供する。

対応しているソルバーは以下の通り。

  • gurobi

  • cplex

  • copt

  • xpress

  • mosek

  • highs

  • cbc

  • scip

  • gcg

MPインターフェイスによって、ソルバーが対応していない場合には、以下を自動的に再定式化する。

  • 条件演算子: if-then-else; ==>, <==, <==>

  • 論理演算子: or, and, not; exists, forall

  • 区分線形関数: abs; min, max; <<breakpoints; slopes>>

  • 計数演算子: count; atmost, atleast, exactly; numberof

  • 関係演算子および比較演算子: >(=), <(=), (!)=; alldiff

  • 相補性演算子: complements

  • 非線形演算子および関数: *, /, ^; exp, log; sin, cos, tan; sinh, cosh, tanh

  • 集合所属演算子: in

離接制約 or

reset;
var x >= -1000 <= 1000;
var y >= -1000 <= 1000;
maximize total: 5 * x + 2 * y;
s.t. only_one_positive: x <= 0 or y <= 0;

option solver highs;
solve;
display x,y;

solexpand only_one_positive;
HiGHS 1.10.0: HiGHS 1.10.0: optimal solution; objective 5000
0 simplex iterations
0 branching nodes
x = 1000
y = 0

subject to only_one_positive:
    x <= 0 || y <= 0;

if then ==>

reset;
var x >= -1000 <= 1000;
var y >= -1000 <= 1000;
maximize total: 5 * x + 2 * y;
s.t. only_one_positive: x > 0 ==> y <= 0;

option solver highs;
solve;
display x,y;

solexpand only_one_positive;
HiGHS 1.10.0: HiGHS 1.10.0: optimal solution; objective 5000
0 simplex iterations
0 branching nodes
x = 1000
y = 0

subject to only_one_positive:
    !(x > 0) || y <= 0;

if and only if <==>

reset;
var x >= -1000 <= 1000;
var y >= -1000 <= 1000;
minimize total: 5 * x + 2 * y;
s.t. exactly_one_positive: x > 0 <==> y <= 0;

option solver highs;
solve;
display x,y;
HiGHS 1.10.0: HiGHS 1.10.0: optimal solution; objective -4999.9998
0 simplex iterations
0 branching nodes
x = -1000
y = 0.0001

選言標準形 disjunctive normal form

reset;

var x >= -1000 <= 1000;
var y >= -1000 <= 1000;
minimize total: 5 * x + 2 * y;
s.t. exactly_one_positive_with_gap:
    (x <= 0 and y >= 3) or (x >= 3 and y <= 0);

option solver highs;
solve;
display x,y;

solexpand exactly_one_positive_with_gap;
HiGHS 1.10.0: HiGHS 1.10.0: optimal solution; objective -4994
0 simplex iterations
0 branching nodes
x = -1000
y = 3

subject to exactly_one_positive_with_gap:
    x <= 0 && y >= 3 || x >= 3 && y <= 0;

全相異制約 alldiff

8-クイーンの例

reset;
param n integer > 0;             # N queens
var Row {1..n} integer >= 1 <= n;
s.t. row_attacks: alldiff ({j in 1..n} Row[j]);
s.t. diag_attacks: alldiff ({j in 1..n} Row[j]+j);
s.t. rdiag_attacks: alldiff ({j in 1..n} Row[j]-j);

let n := 8;
option solver gurobi;
solve;
display Row;

option gurobi_auxfiles rc; # export row/col names
option gurobi_options 'writegraph=model.jsonl';
Gurobi 12.0.1: Gurobi 12.0.1: optimal solution
0 simplex iterations
Objective = find a feasible point.
Row [*] :=
1  5
2  2
3  4
4  7
5  3
6  8
7  6
8  1
;

if .. then ..

施設配置問題の例

reset;
# Set up the sets
set I := 1..5;  # potential facility locations
set J := 1..10;  # customers

# Set up the parameters
param w {i in I} = Normal(60, 20);  # fixed costs for each facility
param u {i in I, j in J} = Uniform(10, 30);  # transportation costs from each facility to each customer
param d {j in J} = Uniform(5, 10);  # demand for each customer

# Set up the decision variables
var x {i in I, j in J} >= 0 <= d[j];  # amount of demand for customer j satisfied by facility i

# Set up the objective function
minimize total_cost:
    sum {i in I, j in J} u[i,j]*x[i,j] +
        sum {i in I} if sum {j in J} x[i,j] > 0 then w[i];

# Set up the constraints
subject to demand_constraint {j in J}:
    sum {i in I} x[i,j] = d[j];

option solver highs;
solve;
display x;
HiGHS 1.10.0: HiGHS 1.10.0: optimal solution; objective 1216.446474
7 simplex iterations
1 branching nodes
x [*,*] (tr)
:       1         2      3      4         5       :=
1    0         0         0   0         7.00354
2    0         7.6232    0   0         0
3    9.89917   0         0   0         0
4    0         0         0   0         9.71021
5    0         0         0   6.95342   0
6    0         5.3371    0   0         0
7    0         0         0   0         5.71001
8    0         0         0   8.05193   0
9    9.23089   0         0   0         0
10   0         9.11613   0   0         0
;

不連続な変数の領域

栄養問題の例

reset;
model;
set NUTR;
set FOOD;
param cost {FOOD} > 0;
param f_min {FOOD} >= 0;
param f_max {j in FOOD} >= f_min[j];
param n_min {NUTR} >= 0;
param n_max {i in NUTR} >= n_min[i];
param amt {NUTR,FOOD} >= 0;

# 4. A union of two intervals
var Buy {FOOD} in interval[2, 5] union interval[7,10];
#var Buy {FOOD} in {0} union interval [10, 50];#
#var Buy {FOOD} in 0..4 union 9..13 union 17..20;

minimize Total_Cost: sum {j in FOOD} cost[j] * Buy[j];
subject to Diet {i in NUTR}:
n_min[i] <= sum {j in FOOD} amt[i,j] * Buy[j] <= n_max[i];

data;
set NUTR := A B1 B2 C ;
set FOOD := BEEF CHK FISH HAM MCH MTL SPG TUR ;
param:   cost  f_min  f_max :=
BEEF   3.19    0     100
CHK    2.59    0     100
FISH   2.29    0     100
HAM    2.89    0     100
MCH    1.89    0     100
MTL    1.99    0     100
SPG    1.99    0     100
TUR    2.49    0     100 ;
param:   n_min  n_max :=
A      700   10000
C      700   10000
B1     700   10000
B2     700   10000 ;
param amt (tr):
        A    C   B1   B2 :=
BEEF   60   20   10   15
CHK     8    0   20   20
FISH    8   10   15   10
HAM    40   40   35   10
MCH    15   35   15   15
MTL    70   30   15   15
SPG    25   50   25   15
TUR    60   20   15   10 ;

option solver highs;
solve;
display Buy;
HiGHS 1.10.0: HiGHS 1.10.0: optimal solution; objective 101.0133333
7 simplex iterations
1 branching nodes
Buy [*] :=
BEEF   2
 CHK  10
FISH   2
 HAM   2
 MCH  10
 MTL  10
 SPG   7.33333
 TUR   2
;

区分線形関数

### 絶対値 abs

reset;

var x {1..2} >=-30 <=100;
minimize Objective: abs(x[1]) - 2*abs(x[2]);
s.t. Con1: 3*x[1] - 2*x[2] <=  8;
s.t. Con2: x[1] +   x[2] == 14;

option solver highs;
solve;
display x;

#solexpand Objective;
HiGHS 1.10.0: HiGHS 1.10.0: optimal solution; objective -58
1 simplex iterations
1 branching nodes
x [*] :=
1  -30
2   44
;

min, max

reset;

var x {1..2} >=0;
minimize Objective: max( 3 * x[1] + 4 * x[2], 2 * x[1] + 7 * x[2]);
s.t. Con1: x[1] + 2 * x[2] >= 12;
s.t. Con2: 2 * x[1] + x[2] >= 15;

option solver highs;
solve;
display x;
HiGHS 1.10.0: HiGHS 1.10.0: optimal solution; objective 31.2
4 simplex iterations
0 barrier iterations
x [*] :=
1  7.2
2  2.4
;

reset;

var x {1..2} >=0, <=1000; #Big Mを使うので、明示的に境界を設定する； gurobiだと設定なしでも解けるが、heighsではエラーする、
minimize Objective: min( 3 * x[1] + 4 * x[2], 2 * x[1] + 7 * x[2]);
s.t. Con1: x[1] + 2 * x[2] >= 12;
s.t. Con2: 2 * x[1] + x[2] >= 15;

option solver highs;
solve;
display x;
HiGHS 1.10.0: HiGHS 1.10.0: optimal solution; objective 24
8 simplex iterations
1 branching nodes
x [*] :=
1  12
2   0
;

区分線形関数

def f(x):
  return x**3-2*x**2+5

def fprime(x):
  return 3*x**2-4*x

xlist=[0,1,2,3]
ylist=list(map(f,xlist))
slope = [fprime(xlist[0])]
for i,x in enumerate(xlist[:-1]):
  slope.append(f(xlist[i+1])-f(x))
slope.append(fprime(xlist[-1]))
print(slope)
[0, -1, 1, 9, 15]
reset;
var x >=0, <=1000;
param xlist{0..3};
param slope{0..4};
minimize Objective: << {i in 0..3} xlist[i];  {i in 0..4} slope[i] >> x ;
ampl.param["xlist"] = xlist
ampl.param["slope"] = slope
display slope;
option solver highs;
solve;
display x;
slope [*] :=
0   0
1  -1
2   1
3   9
4  15
;

HiGHS 1.10.0: HiGHS 1.10.0: optimal solution; objective -1
0 simplex iterations
0 barrier iterations
x = 1

数え上げ count, atleast, atmost

輸送問題を例とする。

reset;
model;
set ORIG;   # origins
set DEST;   # destinations

param supply {ORIG} >= 0;   # amounts available at origins
param demand {DEST} >= 0;   # amounts required at destinations

   check: sum {i in ORIG} supply[i] = sum {j in DEST} demand[j];

param cost {ORIG,DEST} >= 0;   # shipment costs per unit
var Trans {ORIG,DEST} >= 0;    # units to be shipped

minimize Total_Cost:
   sum {i in ORIG, j in DEST} cost[i,j] * Trans[i,j];

subject to Supply {i in ORIG}:
   sum {j in DEST} Trans[i,j] = supply[i];

subject to Demand {j in DEST}:
   sum {i in ORIG} Trans[i,j] = demand[j];

#少なくとも7本の経路が600以上の輸送量がある
subject to Min_Serve:
     count {i in ORIG, j in DEST} (Trans[i,j] >= 600) >= 7;

#高々2本の経路が1000以上の輸送量がある（等号のtrickで3本のように見えるが、1000-epsilonと解釈）
#subject to at_most_k:
#    atmost 2 {i in ORIG, j in DEST} (Trans[i,j] >= 1000);

data;
param: ORIG:  supply :=  # defines set "ORIG" and param "supply"
        GARY   1400
        CLEV   2600
        PITT   2900 ;

param: DEST:  demand :=  # defines "DEST" and "demand"
        FRA     900
        DET    1200
        LAN     600
        WIN     400
        STL    1700
        FRE    1100
        LAF    1000 ;

param cost:
         FRA  DET  LAN  WIN  STL  FRE  LAF :=
   GARY   39   14   11   14   16   82    8
   CLEV   27    9   12    9   26   95   17
   PITT   24   14   17   13   28   99   20 ;

option solver gurobi; #highsだと解けない
solve;
display Trans;
Gurobi 12.0.1:   tech:writegraph = model.jsonl
Gurobi 12.0.1: optimal solution; objective 197000
26 simplex iterations
1 branching node
Trans [*,*] (tr)
:     CLEV   GARY   PITT    :=
DET   1200      0      0
FRA      0      0    900
FRE      0   1100      0
LAF      0    300    700
LAN    600      0      0
STL    600      0   1100
WIN    200      0    200
;

非線形関数

  • log (expr), log10 (expr)
  • exp (expr)
  • sin (expr), cos (expr), tan (expr), asin (expr), acos (expr), atan (expr)
  • sinh (expr), cosh (expr), tanh (expr), asinh (expr), acosh (expr), atanh (expr)
  • expr1 ^ expr2

ソルバーが対応可能な場合はそのまま、非対応の場合には区分的線形関数に近似する。

例として、3次関数 \(x^3+2x^2+5\) の最小化を行う。区分的線形近似のためには、変数 \(x\) の範囲を明示的に定義しておく必要がある。

reset;
model;
var x >=0, <=3;
minimize Objective: x^3-2*x^2+5;
option solver scip; #highsだと間違える；gurobiかscipもしくは非線形最適化ソルバーを用いる
solve;
display x;
SCIP 9.0.1: WARNING: No dual information available when presolving was performed.
SCIP 9.0.1: optimal solution; objective 3.814813982
169 simplex iterations
x = 1.33352