この記事の内容

ソルバーしか使えない私が、栄養最適化問題で、食品の種類を最小にしたくて、非線形数理最適化をやってみた結果を紹介します。

要素数最小化の意義

例えば栄養最適化の場合、20種類もの食品の組み合わせになると買うのが面倒です。
個人の買い物くらいならまだいいですが、これが栄養最適化食品工場のだと考えると、原料の仕入れ、在庫の調整や生産に使うラインの数を食品ごとに用意することの運用コストの方が原料コストを上回る気がします。

最適化問題

通常の栄養最適化問題は、必要な栄養摂取量の範囲に収まる中で最安値となる食品の組み合わせを探す問題になります。
これは整数線形計画の問題として書けて、
制約条件:
Min(N) < sum( N(f) * x(f) ) < Max(N)
目的関数:
sum( cost(f) * x(f) )

N(f): 栄養素 N について、食品 f に含まれる栄養素の量
x(f): 食品 f の個数
Min(N), Max(N): 栄養素 N の摂取量の上下限値
cost(f)：食品 f の1個当たりの値段

食品の「種類数」を最小化しようとすると、食品 f を使うかどうかの変数a(f)∈(0,1) を使って
sum(a(f))を最小化すればよいのですが、その際の制約条件が
Min(N) < sum( N(f) * a(f) * x(f) ) < Max(N)
となり、線形式で表すことができなくなります。

食品の種類を要素数の最小化は線形計画法では解けないことが分かったので、非線形整数計画を使うことにします。

制約条件は同じで、目的関数を変えます。目的関数はx(f)が、0なら0になり、1以上なら1に近くなる関数で合計値最小化にする方式で考えました。
個数だけなのでシンプルに
目的関数：
sum(1-100^(-x(f))
とします。

なお、線形計画法でやるなら、用意したすべての食品について、2個、3個と、数の少ない順でそれぞれの組み合わせについてそれぞれ線形計画問題を解き、それで解が出た個数が答えになります。(非線形計画法を使うよりそっちの方が早いと思っていますがまだ試したことはありません)

環境構築(Windows)

Windows10
Python 3.8.1 64bit
Bonmin

フリーの非線形整数計画(MINLP)のソルバーとして、Coin-ORプロジェクトで開発された凸最適化MINLPソルバーのBonminを使います。
凸最適化についてはWikipediaなどに詳しいですが、上記の制約条件、目的関数ともに凸集合、凸関数なので凸最適化ソルバーのBonminを使うことができます。
Bonminは残念ながら最適な解(グローバル解)は保証されていません。(同じくCoin-ORプロジェクトのCouenneは非凸最適化問題のソルバーですが、こちらはグローバル解を求められます)

WindowsでPythonとBonminを使って環境構築

AMPLからダウンロード

最適化モデリング言語のAMPLのオープンソースソルバーのダウンロードページから最新のWindows64bit版のBonminをダウンロードできます（2022年8月時点で1.8.8）
https://ampl.com/products/solvers/open-source/
COIN-OR公式Githubで紹介されているので安心して使えます。

※COIN-ORのプロジェクトページからもダウンロードできますが、古いです。(https://www.coin-or.org/download/binary/Bonmin/ からBonmin-1.4.0-win32-msvc9.zip をダウンロード可能)

Zipフォルダを展開

AMPLからダウンロードしたZipファイルを展開します。
bonmin-win64\bonmin.exe
bonmin-win64\coin-license.txt
bonmin-win64\libipoptfort.dll

Python実行環境の準備

bonmin-win64と同じフォルダにpythonの環境bonminを作成し、必要なライブラリをインストール
python -m venv bonmin
./bonmin/Scripts/Activate
pip install pyomo numpy

Dockerで環境構築(Windows,Linux,MacOS）

下記を参考にしてみてください
https://hub.docker.com/r/coinor/coin-or-optimization-suite

コード

bonmin-win64と同じフォルダに下記のbonmintest.pyを保存

bonmintest.py

from pyomo.environ import *
from pyomo.opt import SolverFactory
import numpy as np

F = 5 #Foodの数
N = 3 #栄養素の数
food_data = np.zeros((F,N))

# 食品に含まれる栄養素
food_data[0] = [1, 0, 0] #food1の100gあたりの各栄養素
food_data[1] = [0, 1, 0] #food2の100gあたりの各栄養素
food_data[2] = [0, 0, 1] #food3の100gあたりの各栄養素
food_data[3] = [3, 6, 1] #food4の100gあたりの各栄養素
food_data[4] = [1, 2, 1] #food5の100gあたりの各栄養素

# 摂取量の上下限
nutrient_max = [10, 20, 5] # 各栄養素の一日の最大値
nutrient_min = [3, 6, 2] # 各栄養素の一日の最小値


M = ConcreteModel("Bonmin Test") #モデル
M.foods = Set(initialize=range(F)) #変数の数
M.var = Var(M.foods, within=NonNegativeIntegers) #最適化対象の変数定義：非負整数NonNegativeIntergersを指定。

# 制約条件を返す関数
def _con1(model,foods,nutrient):
    return (nutrient_min[nutrient],sum(food_data[food][nutrient]*model.var[food] for food in foods),nutrient_max[nutrient])

# 制約条件
M.constraints = ConstraintList()
for nutrient in range(N):
    M.constraints.add(_con1(M,M.foods,nutrient))

# 目的関数を返す関数
def obj_rule(model):
    return sum((1-100**(-model.var[food])) for food in model.foods)

# 最小化問題として目的関数を設定
M.value = Objective(rule=obj_rule, sense=minimize)

# ソルバーとしてBonminを指定
opt = SolverFactory("bonmin", executable="./bonmin-win64_1.8.8fromAMPLcom/bonmin")

# オプションの設定 https://www.coin-or.org/Bonmin/option_pages/options_list_bonmin.html
bonmin_option ={"bonmin.algorithm":"B-BB",#B-BB,B-OA,B-QG,B-Hyb,B-Ecp,B-iFP,Cbc_Par
                #"bonmin.pump_for_minlp":"yes", # whether to run the feasibility pump heuristic for MINLP
                "bonmin.random_generator_seed":-1,#-1 sets seeds to time since Epoch/ option for B-BB
                #"bonmin.time_limit":10,
                #"bonmin.cutoff":8,
                "halt_on_ampl_error":"yes"}

# 最適化モデルの確認
print(M.constraints.pprint())

# 実行
result = opt.solve(M, options = bonmin_option, tee=True)

# 結果の表示
M.display()

# 結果の確認
condition = result["Solver"][0]['Termination condition']
print(condition) #optimal=最適化できた。infeasible=解無し

# 結果の抽出
if condition == "optimal":
    print("選択した食品の種類:{}種類".format(round(M.value())))

    print("選択した食品と個数")
    for f,food in zip(M.var,range(F)):
        var = M.var[f]()
        print("food{}:{}個".format(f,var))

    print("取得できる栄養素（下限～上限）")
    for const in M.component_objects(Constraint):
        for c in const:
            print("栄養素{}={}({}~{})".
                  format(c,M.constraints[c](),
                         value(M.constraints[c].lower),
                         value(M.constraints[c].upper)))

これを
python bonmintest.py
として実行すると下記の結果が得られると思います

constraints : Size=3, Index=constraints_index, Active=True
    Key : Lower : Body                             : Upper : Active
      1 :   3.0 :     var[0] + 3.0*var[3] + var[4] :  10.0 :   True
      2 :   6.0 : var[1] + 6.0*var[3] + 2.0*var[4] :  20.0 :   True
      3 :   2.0 :         var[2] + var[3] + var[4] :   5.0 :   True
None
Bonmin 1.8.8 using Cbc 2.10.5 and Ipopt 3.12.13
bonmin: bonmin.algorithm=B-BB
bonmin.random_generator_seed=-1
halt_on_ampl_error=yes
bonmin.algorithm=B-BB
bonmin.random_generator_seed=-1
halt_on_ampl_error=yes


******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

NLP0012I 
              Num      Status      Obj             It       time                 Location
NLP0014I             1         OPT 0.99999981       55 0.06
NLP0012I 
              Num      Status      Obj             It       time                 Location
NLP0014I             1         OPT 0.99999882       19 0.023
NLP0014I             2         OPT 0.99999882        9 0.009
NLP0014I             3         OPT 0.999999        0 0
NLP0012I
              Num      Status      Obj             It       time                 Location
NLP0014I             1         OPT 0.999999        0 0
Cbc0012I Integer solution of 0.999999 found by DiveMIPFractional after 0 iterations and 0 nodes (0.03 seconds)      
Cbc0001I Search completed - best objective 0.999999, took 0 iterations and 0 nodes (0.03 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost

        "Finished"
Model Bonmin Test

  Variables:
    var : Size=5, Index=foods
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :   0.0 :  None : False : False : NonNegativeIntegers
          1 :     0 :   0.0 :  None : False : False : NonNegativeIntegers
          2 :     0 :   0.0 :  None : False : False : NonNegativeIntegers
          3 :     0 :   0.0 :  None : False : False : NonNegativeIntegers
          4 :     0 :   3.0 :  None : False : False : NonNegativeIntegers

  Objectives:
    value : Size=1, Index=None, Active=True
        Key  : Active : Value
        None :   True : 0.999999

  Constraints:
    constraints : Size=3
        Key : Lower : Body : Upper
          1 :   3.0 :  3.0 :  10.0
          2 :   6.0 :  6.0 :  20.0
          3 :   2.0 :  3.0 :   5.0
optimal
選択した食品の種類:1種類
選択した食品と個数
food0:0.0個
food1:0.0個
food2:0.0個
food3:0.0個
food4:3.0個
取得できる栄養素（下限～上限）
栄養素1=3.0(3.0~10.0)
栄養素2=6.0(6.0~20.0)
栄養素3=3.0(2.0~5.0)

まとめ

Bonminを用いて非線形整数計画（MINLP）を解くことができました。
この例では、最初にも書いたように、要素の組み合わせごとに線形計画問題を解いていけば同じ結果が得られます（food4一つで制約を満たせる、という結果がすぐ出ます）

それでは、非線形で何が嬉しいかというと、数式一つで「要素数とコストのバランスがよい解を探す」といったことが可能になります。要素は減らせたが原料が高い場合、もう一つ増やしてでも安くできた方が嬉しい（ライン増設のコストと原料コストを比較して最適化）わけです。（規模が大きくなると時間がかかるので、順次線形計画法で計算した方が早いかもしれませんが…）

Bonminを使ってもっと大規模な栄養最適化問題を解くアプリを作ってみたので、よかったら見てくださいDietaryOptimizer。

参考サイト

凸最適化(Wikipedia)
凸最適化(Convex Optimization)の基礎 - MyEnigma (hatenablog.com)
Convex Optimization(PDF)

COIN-OR(Github.com)
COIN-OR(Github.io)
Bonmin
Bonmin Options List
CasADiとBONMINで混合整数非線形計画問題を解く(Helveさん)

Pyomo
Pyomo(制約条件の取得)