MATLAB Optimization ToolboxでMILPやMIQPの最適化問題を定式化してCPLEX v20.1で解く

普段からよく混合整数計画問題（MILP）や混合整数二次計画問題（MIQP）の最適化を取り扱っている者です。MATLAB Optimization Toolbox ¹の問題ベースによる最適化問題の定式化は設定がしやすく，ライブスクリプトで各制約の上に式を書くことができるので非常に取り扱いやすいので多用しています。

しかしながら，問題の規模にもよりますが，MATLABの最適化ソルバー(intlinprog,など）はIBM社のCPLEX ² のようなサードパーティー製の有償ソルバーには及ばず，計算時間がかなり違ったり，解が出てこないこともあります。

CPLEXではこれまでMATLAB Toolbox APIがありましたので，CPLEXをインストールすれば，MATLABの関数として呼び出すことができましたが，2020年の最新バージョン(v20.1)では，MATLAB Toolbox APIが削除されましたので，MATLABのコマンドからはそのまま使うことできなくなりました。そのため，

• 古いバージョンを使い続ける
• 最適化計算は別のプログラム環境で実施する

といったことを考えましたが，CPLEXの最新バージョンの計算速度の向上が大きかったので，後者を実施することにしました。また，MATLABでの実行を前提としていますので，得られた解もMATLABの変数に入れることを想定しています。

仕様

今回は別のプログラム環境として，CPLEXがサポートしているPython 3.7を用いました。
なお，CPLEXおよびpythonの使用環境は構築済み（それぞれが同一フォルダから呼び出し可）であることを前提としています。

• MATLAB 2020b (Optimization Toolbox)
• CPLEX v12.10以下（MATLABのフォルダパスに追加しておきます）
• CPLEX v20.1（新バージョンで解く場合）
• python 3.7.9

MATLABコード

**prob2cplexOpt.m**

prob2cplexOpt.m

strProb.Options = cplexoptimset('timelimit',1);     %最大計算時間
strProb.Options.exportmodel = 'myModel.lp';         %lpファイル名
strProb.ctype=[];                                   %CPLEX用のx変数の型の設定
strProb.ctype=repelem('C',length(strProb.ub));      %CPLEX用のx変数の型 C:連続値
strProb.ctype(1,strProb.intcon(1,:))='I';           %CPLEX用のx変数の型 I:整数値
[x] = cplexmiqp (strProb);                          %CPLEXの実行とlpファイルの出力
!python cplexOpt.py                                 %pythonの実行   
x=csvread('x_sol.csv',1,1);                         %csvファイルの読み込み

Pythonコード

**cplexOpt.py**

cplexOpt.py

from __future__ import print_function
import cplex
import csv  
import array
import pandas as pd

def optmization():
    
    prob = cplex.Cplex()     
    prob.read("myModel.lp") #lpファイルの読み込み 
    prob.solve()
    return prob

def write(prob):
    sol = prob.solution

    # solution.get_status() returns an integer code
    print("Solution status = ", sol.get_status(), ":", end=' ')
    # the following line prints the corresponding string
    print(sol.status[sol.get_status()])
    print("Solution value  = ", sol.get_objective_value())

    xList =[]
    x0=0
    numcols = len(sol.get_values())
    for j in range(numcols):
             xNo='x'+str(j+1)
             x0=sol.get_values(xNo)
             xList.append([x0])

    # Pandasのデータフレームにリストを格納
    df = pd.DataFrame(xList,columns=['x'])
    # 書き出し
    df.to_csv('x_sol.csv', mode='w')

if __name__ == "__main__":
    prob=optmization()
    write(prob)

解説・手順

以下の手順でプログラムを実行します。

1. MATLAB のOptimization Toolboxで解きたい最適化問題を定式化
2. Optimization Toolbox内の関数 prob2structによりソルバーベースの問題構造体へ変換
3. CPLEX特有の必要なオプションの設定追加
4. MATLAB上でCPLEX v12.10 で実行し，lp(あるいはmps)ファイル生成
5. Python上でCPLEX v20.1を実行

手順1:MATLAB で最適化問題を定式化

ここでは省略しますが，問題ベースの最適化の設定を参考に最適化問題を定式化します。
※希望があれば、ライブスクリプトに問題ベースの最適化を行った例を公開します。

手順2:ソルバーベースの問題構造体へ変換

Optimization Toolbox内の関数prob2structによりソルバーベースの問題構造体strProbへ変換します。
「なぜわざわざソルバーベースから始めないのか」といった声が聞こえてきそうですが、視覚的にわかりやすいから、手順1のライブスクリプトに問題ベースを定式化します。

prob2cplexOpt.m

strProb=prob2struct(problem)

手順3:CPLEX特有の必要なオプションの設定追加

ソルバーベースしたstrProbにCPLEX特有のオプションを追加します。
cplexoptimsetでCPLEXのオプションを取得しますが，MATLABのCPLEX v12.10の実行ではlpファイルを生成することが目的ですので，計算時間は1秒にして，即計算が終わるようにします。また，lpファイル名を指定しておきますが，拡張子を.mpsにするとmpsファイルを自動的に生成してくれます。

prob2cplexOpt.m

strProb.Options = cplexoptimset('timelimit',1);     %最大計算時間
strProb.Options.exportmodel = 'myModel.lp';         %lpファイル名

また，CPLEXではctypeで決定変数の型（整数型・連続型等）を文字で設定する必要があります。ここでは，大多数が連続値（'C'）であり，一部が整数であるとし，整数の要素を示すstrProb.intconに従って整数型('I')を代入させていきます。

prob2cplexOpt.m

strProb.ctype=[];                                   %CPLEX用のx変数の型の設定
strProb.ctype=repelem('C',length(strProb.ub));      %CPLEX用のx変数の型 C:連続値
strProb.ctype(1,strProb.intcon(1,:))='I';           %CPLEX用のx変数の型 I:整数値

手順4:MATLAB上でCPLEX v12.10 で実行し，lp(あるいはmps)ファイル生成

MATLAB上でCPLEXを実行してlpファイルを生成します。ちなみに，MILPでも二次目的関数

min \bigl( \frac{1}{2} \boldsymbol{x}'\boldsymbol{H}\boldsymbol{x}+\boldsymbol{f}\boldsymbol{x}  \bigr) \\
ただし，\boldsymbol{x} =\bigl(x_{1},x_{2},...,x_{n}\bigr)

で$\boldsymbol{H}=0$ として取り扱えるので，私はcplexmiqpを使っています。

prob2cplexOpt.m

[x] = cplexmiqp (strProb);                          

なお，ここでは最適化時間は1秒としていますので，解が得られないことが多いでしょう。逆に，それくらいの時間で解けるのであれば，ここで得られるxが解として終了して問題ないと思います。解けない場合は次の手順に進みましょう。

手順5:Python上でCPLEX v20.1を実行

pythonによりCPLEX v20.1による最適化を実行をします。
pythonのコードは以下のようにして構築しますが，インポートはcplex，csv，array，pandasなどをしておきます。
ここではざっくりと説明しますが，

• 関数optmizationでMATLABで生成したlpファイルを読み込み，最適化
• 関数writeで最適化結果をcsvファイルに出力

しています。

なお，lpファイル内の変数（$x1,x2,x3,...$）の番号と得られる解の順番は対応していません。
これはlpファイルから読み取る変数（$x1,x2,x3,...$）は，全て単なる文字列として扱われるためであり，lpファイルの上から勝手に番号が振られます。

Minimize
 obj1: x4 + x3 + x2 + x1

と書いてある場合はCPLEXの解$\boldsymbol{x}$は（$x4,x3,x2,x1$）の順になります。そうすると，MATLABの変数'x'とインデックスが異なってしまうので，xの変数文字列xNoで$x1,x2,x3,x4,...$となるようにソートします。

cplexOpt.write

    xList =[]
    x0=0
    numcols = len(sol.get_values())
    for j in range(numcols):
             xNo='x'+str(j+1)
             x0=sol.get_values(xNo)
             xList.append([x0])

ソートしたら，そのリストをcsvファイル出力します。解以外にも最適化のステータスも出力しても良いとは思いますが，今回のpythonのコードはここまでにしています。
後はMATLABあるいはpythonで実行するだけです。MATLABで実行する場合は以下のように!で始めれば実行でき，csvファイルを読み込みすればよいです。

prob2cplexOpt.m

!python cplexOpt.py                                 %pythonの実行   
x=csvread('x_sol.csv',1,1);                         %csvファイルの読み込み

計算時間・補足等

決定変数が少ない場合は，MATLABのintlinprogで解けますし，CPLEX v12.10でも十分なのでここまでする必要はないかと思います。私が現在扱っている決定変数が21万個，制約条件も35万個くらいある大規模なMILPですと，手順1～5の全ての手順の時間を含めても計算時間は3～４割くらい短くなりました。参考になれば幸いです。

手法	計算時間
intlinprog	580秒
CPLEX v12.10(MATLABのみ)	240秒
CPLEX v20.1 (MATLAB→python)	160秒

1. MATLABのOptimization Toolboxの詳細はこちら ↩

2. CPLEXの詳細はこちらから見れます。また，CPLEXのアカデミック版は無料で使用することができます。 ↩