見出し画像

ExcelとLpSolve

Excelと斬鉄剣

目的

最適化問題の解法としてよく用いられる混合整数計画法（MIP）。MIPはソルバーでも解くことができますが、ソルバーには決定変数の数が200まで（200変数を超えるとエラーになる）という制約があるので、無償で利用できるLpSolveをExcelから利用できるようにしておくとなにかと便利です。

要件

ユーザインタフェース

demo.xlsのサンプルを少し変更して次のような混合整数計画問題を考えます。
決定変数：x1 , x2 , x3 , x4 , x5（非負）
目的関数：2 x1 + 3 x2 - 2 x3 + 3 x4 + 3 x5 → 最大化
制約条件
　　　c1：0.5 x1 + 4 x2 + x3 + 2 x4 + 3 x5 >= 3
　　　c2：x1 + 2 x2 + x3 + x4 + 2 x5 = 8
　　　c3：4 x1 + x2 - 2 x3 + 2 x4 - 2 x5 <= 5
　　　x4は整数、x5はバイナリ
　　　x5 <= 5.3 , x2 >= 2
この問題を記述できるユーザインタフェースとして次のようなシートを用意します。水色ハッチング部分は変数名および制約条件名ですがおそらく日本語に対応していないので半角英数で記入します。黄色ハッチング部分には混合整数計画問題を解いたあとに最適化結果である決定変数や目的関数の値を入力します。

実行前
実行後

決定変数の列や制約条件の行を挿入すればある程度問題の拡張にも対応可能です。

別の問題の実行後

プログラムソース

demo.xlsのサンプルソースはプログラム内に問題をべた書きされているので、変数の数や制約条件の数、最大化か最小化か、変数の上下限、変数のタイプ（整数やバイナリ）は画面から読み取るように組み込みます。cは変数の数のカウンタ、c1は制約の数のカウンタです。制約条件の結果も取得していますが、画面に書き込んではいません。

Option Explicit

Sub 混合整数計画法を解く()

    Dim lpsolve As lpsolve55
    
    Const 最大行 = 1048576
    Const 最大列 = 16384
    
    Dim i As Long
    Dim j As Long
    Dim c As Long
    Dim c1 As Long
    Dim t As String
    Dim t1 As String
    Dim a() As Double
    Dim lp As Long
    Dim obj As Double
    Dim st As Long
    Dim row() As Double
    Dim col() As Double
    Dim 目的行 As Long
    Dim 結果行 As Long
    Dim ws As Worksheet
    Set ws = Sheets("問題")
    ws.Activate
    
    '変数のカウント------------------------------------------------------------
    c = 0
    For i = 2 To 最大列
        If ws.Cells(1, i).value = "符号" Then
            Exit For
        Else
            c = c + 1
        End If
    Next i
    
    'lpsolveの初期化-----------------------------------------------------------
    Set lpsolve = New lpsolve55
    lpsolve.Init Application.ActiveWorkbook.Path
    lp = lpsolve.make_lp(0, c)
    
    '画面設定の読込------------------------------------------------------------
    
    c1 = 0
    For i = 1 To 最大行
        t = ws.Cells(i, 1).value
        t1 = ws.Cells(i, c + 2).value
        
        If t = "" Then
            Exit For
        
        '変数名の設定
        ElseIf t = "変数" Then
            For j = 1 To c
                lpsolve.set_col_name lp, j, ws.Cells(i, j + 1).value
            Next j
            
        '目的関数の設定
        ElseIf t = "目的関数" Then
            目的行 = i
            ReDim a(0 To c) As Double
            For j = 1 To c
                a(j) = Val(ws.Cells(i, j + 1).value)
            Next j
            lpsolve.set_obj_fn lp, a(0)
            '最大化する場合
            If t1 = "max" Then
                lpsolve.set_maxim lp
            End If
        
        '制約条件の設定（制約条件名の設定）
        ElseIf InStr(t1, "=") > 0 Then
            c1 = c1 + 1
            ReDim a(0 To c) As Double
            For j = 1 To c
                a(j) = Val(ws.Cells(i, j + 1).value)
            Next j
            If t1 = ">=" Then
                lpsolve.add_constraint lp, a(0), GE, Val(ws.Cells(i, c + 3).value)
            ElseIf t1 = "=" Then
                lpsolve.add_constraint lp, a(0), EQ, Val(ws.Cells(i, c + 3).value)
            ElseIf t1 = "<=" Then
                lpsolve.add_constraint lp, a(0), LE, Val(ws.Cells(i, c + 3).value)
            End If
            lpsolve.set_row_name lp, c1, ws.Cells(i, 1).value
        
        'タイプの設定
        ElseIf t = "タイプ" Then
            For j = 1 To c
                If ws.Cells(i, j + 1).value = "int" Then
                    lpsolve.set_int lp, j, True
                ElseIf ws.Cells(i, j + 1).value = "bin" Then
                    lpsolve.set_binary lp, j, True
                End If
            Next j
            
        '変数の上限設定
        ElseIf t = "上限" Then
            For j = 1 To c
                If Val(ws.Cells(i, j + 1).value) > 0 Then
                    lpsolve.set_upbo lp, j, Val(ws.Cells(i, j + 1).value)
                End If
            Next j

        '変数の下限設定
        ElseIf t = "下限" Then
            For j = 1 To c
                If Val(ws.Cells(i, j + 1).value) > 0 Then
                    lpsolve.set_lowbo lp, j, Val(ws.Cells(i, j + 1).value)
                End If
            Next j
        
        '結果行の取得
        ElseIf t = "結果" Then
            結果行 = i
        End If
    
    Next i
        
    'タイムアウト設定
    lpsolve.set_timeout lp, 0
    
    'デバッグモードON
    lpsolve.set_debug lp, True
        
    'トレースモードOFF
    lpsolve.set_trace lp, False
    
    '結果ファイルの設定
    lpsolve.set_outputfile lp, Application.ActiveWorkbook.Path & "\result.txt"

    '結果ファイルへの問題の書き込み
    lpsolve.print_lp lp
    
    '最適化実行
    st = lpsolve.solve(lp)
    
    '目的関数の結果取得
    obj = lpsolve.get_objective(lp)
    
    '変数の結果取得
    ReDim col(1 To lpsolve.get_Ncolumns(lp))
    lpsolve.get_variables lp, col(1)
    
    '最適解が求まったら
    If st = 0 Then
        '変数の結果書き込み
        ReDim a(1 To c) As Double
        ws.Range(Cells(結果行, 2), Cells(結果行, c + 1)).value = col
        '目的関数の結果書き込み
        ws.Cells(目的行, c + 3).value = obj
    Else
        MsgBox ("最適解が求まりませんでした")
    End If
        
    '制約条件の結果取得
    ReDim row(1 To lpsolve.get_Nrows(lp))
    lpsolve.get_constraints lp, row(1)
    
    '結果ファイルへの書き込み
    lpsolve.print_objective lp
    lpsolve.print_solution lp, 1
    lpsolve.print_constraints lp, 1
    
    'lpファイルとmpsファイルの出力
    lpsolve.write_lp lp, Application.ActiveWorkbook.Path & "\lp.lp"
    lpsolve.write_mps lp, Application.ActiveWorkbook.Path & "\lp.mps"

End Sub

lpファイルを出力しておけばデバッグに役立ちます。

lpファイル

最後に

プログラムソースのコピペすら面倒くさいというものぐさな友人のためにExcelファイル添付しておきます。コーヒー1杯分で。笑

ここから先は

101字 / 1ファイル

¥ 100

この記事が気に入ったらチップで応援してみませんか？