MyEnigma

とあるエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

線形計画法(Linear Programming)の概要とPython, Juliaでの解き方

目次

はじめに

先日、最適化問題の一つの方式として二次計画法の概要と、

Pythonでの解き方を紹介しましたが、

myenigma.hatenablog.com

 

今回はもっと多くの分野で使われている

最適化問題の一つである

線形計画法(Linear Programming)の概要と、

同じくいくつかの言語やライブラリでの

解き方を紹介したいと思います。

 

線形計画法の概要

線形計画法は、

下記の式のように、

変数が複数の一次不等式や等式によって制約が与えられている時に、

ある一次のコスト関数を最大・最小化させる方法です。

f:id:meison_amsl:20161025051050p:plain

 

制約条件の不等式は、

非負の変数であるスラック変数という変数を導入することで、

等式制約に変換することができます。

myenigma.hatenablog.com

 

線形計画法を二次元の図に表すと下記のようになります。

Pは制約条件による、状態ベクトルzが取りうる範囲を表しており、

点線はコスト関数の等高線を表しています。

線形計画法の場合、コスト関数の等高線は直線で表されます。

f:id:meison_amsl:20161203063324p:plain

 

この線形計画法は、非常に多くの学問において利用されており、

工学だけでなく、経済学や経営学でも広く利用されています。

 

特にオペレーションズ・リサーチ(OR)と呼ばれる、

様々な計画に対して最も効率的に行う方法を求める学問では、

この線形計画法が非常に広く使われています。

 

この線形計画法を解く上で重要な考え方として、

それぞれの制約条件によって決定する

変数の可能領域を表す多次元凸多面体の頂点にのみ

最適解がある

ということです。

これは、コスト関数や制約上限が一次関数なので、

極値を多面体内部や辺上で取らないためです。

 

この線形計画法を解く方法としては、

上記の考え方を元にした、

シンプレックス法や

(シンプレックス法に関しては、下記を参照下さい)

主双対内点法(Primal-dual interior-point method)が有名です。

myenigma.hatenablog.com

 

Pythonによる線形計画法の解き方

では、実際にPythonを使って、

線形計画法を解いてみたいと思います。

 

cvxoptを使う方法

線形計画法のソルバーとしては、

先日、二次計画法の際に紹介したPython最適化ソルバー

cvxoptに線形計画法を解く機能があるので、

こちらを利用したいと思います。

myenigma.hatenablog.com

 

解く問題としては、

下記の書籍のP171の例題を

実際にPythonを使って解いてみたいと思います。

まずx1,x2,x3,x4の4つの変数があり、

等式制約は下記の通りとします。

f:id:meison_amsl:20161025091927p:plain

f:id:meison_amsl:20161025092025p:plain

また、各変数はすべて0以上という不等式制約もあるとします。

f:id:meison_amsl:20161025123635p:plain

最後に最終的な目的は下記の評価関数を最大化したいと思います。

f:id:meison_amsl:20161025092239p:plain

 

cvxoptでは、

線形計画法を解く場合はsolvers.lpというメソッドを使います。

lpはコスト関数の最小値を求める関数なので、

今回のような最大値を求める問題の場合は

コスト関数の各変数の符号を逆転することで、

最大値を求めることができます。

(下記のコードでは、行列cの部分がマイナスになっている)

また不等式制約についても、

lpの前提条件と問題の条件が逆なので、

係数がマイナスになっていることに注意が必要です。

下記のコードでは、Gとhが不等号制約、

Aとbが等号制約を表しています。

 

"""
Simple sample code to solve linear programming with cvxopt

author Atsushi Sakai
"""

import numpy

import cvxopt
from cvxopt import matrix

c = matrix(numpy.array([-29.0, -45.0, 0.0, 0.0]))
G = matrix(numpy.diag([-1.0, -1.0, -1.0, -1.0]))
h = matrix(numpy.array([0.0, 0.0, 0.0, 0.0]))
A = matrix(numpy.array([[2.0, 8.0, 1.0, 0.0], [4.0, 4.0, 0.0, 1.0]]))
b = matrix(numpy.array([60.0, 60.0]))

print("c:")
print(c)
print("G:")
print(G)
print("h:")
print(h)
print("A:")
print(A)
print("b:")
print(b)

sol = cvxopt.solvers.lp(c, G, h, A, b)

print("x:", sol["x"])
print("obj:", sol["primal objective"])

 

下記が実際に上記のコードを実行した結果になります。

f:id:meison_amsl:20161025092259p:plain

先程の書籍にある通り、

解の答えはx1=10,x2=5,x3=0,x4=0で最大値f=515なので、

fの符号を反転していることを考えると、

ちゃんと正しい値が計算できていることがわかります。

 

上記のコードは下記でも公開しています。

github.com

 

scipyを使う方法

Pythonの科学技術ライブラリとして有名なscipyにも、

内点法による線形計画法を解く関数であるlinprogがあります。

docs.scipy.org

 

下記はこのlinprogを使って、先程のcvxoptで解いた問題と

同じ問題を解くサンプルコードです。

 

from scipy.optimize import linprog
import numpy as np

c = np.array([-29.0, -45.0, 0.0, 0.0])
A_ub = np.diag([-1.0, -1.0, -1.0, -1.0])
b_ub = np.array([0.0, 0.0, 0.0, 0.0])
A_eq = np.array([[2.0, 8.0, 1.0, 0.0], [4.0, 4.0, 0.0, 1.0]])
b_eq = np.array([60.0, 60.0])

result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq)
print(result)

下記が実行結果です。

f:id:meison_amsl:20200405204413p:plain

先程のcvxoptと同じ結果が得られていることがわかります。

個人的には、scipyは様々なPython環境ですでにインストールされていることが多く、

またnumpyのarrayを直接入力できるので、

cvxoptより、より手軽に線形計画法を解くことができると思います。

 

Juliaにおける線形計画法の解き方

続いて、Juliaという比較的新しいプログラム言語で

線形計画法を解いてみたいと思います。

myenigma.hatenablog.com

今回はJuMPというJuliaの最適化ライブラリを使うことにします。

myenigma.hatenablog.com

シンプルな線形計画問題

下記のシンプルな線形計画問題を解いてみたいと思います。

f:id:meison_amsl:20171016022718p:plain

 

上記の問題を解く、Juliaのコードは下記の通りです。

using JuMP
using Clp

m = Model(solver=ClpSolver())

@variable(m, 0<= x1 <=10)
@variable(m, x2 >=0)
@variable(m, x3 >=0)

@objective(m, Max, x1 + 2*x2 + 5*x3)

@constraint(m, -x1 +  x2 + 3*x3 <= -5)
@constraint(m,  x1 + 3*x2 - 7*x3 <= 10)

print(m)
status = solve(m)

println("Optimal Solutions:")
println("x1 = ", getvalue(x1))
println("x2 = ", getvalue(x2))
println("x3 = ", getvalue(x3))

github.com

 

上記を実行すると、下記のように線形計画問題を解くことができます。

f:id:meison_amsl:20171016023033p:plain

 

シンプルな線形計画問題2

以前、cvxpyで解いたシンプルな線形計画問題をJumpで解いてみました。

myenigma.hatenablog.com

   

using JuMP
using Mosek

function main()
    println("JuMP sample2")
    model = Model(solver=MosekSolver())

    @variable(model, 2.5 <= z1 <= 5.0)
    @variable(model, -1.0 <= z2 <= 1.0)
    @NLobjective(model, Min, abs(z1+5.0) + abs(z2-3.0))

    println("The optimization problem to be solved is:")
    println(model) # Shows the model constructed in a human-readable form

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(z2))
end

main()

github.com

 

上記を実行すると下記のように、

cvxpyで同じ問題を解いた時と同じ結果が得られます。

f:id:meison_amsl:20170731012256p:plain

 

一つ注意点として、abs関数をコスト関数の中で利用する場合は、

@objectiveではなく、@NLobjectiveを使う必要があるみたいです。

stackoverflow.com

 

下記はその他のサンプルです。

f:id:meison_amsl:20170801064857p:plain

using JuMP
using Ipopt

function main()
    println("JuMP sample3")
    model = Model(with_optimizer(Ipopt.Optimizer))

    @variable(model, z1>=0.0)
    @variable(model, z2>=0.0)
    @objective(model, Min, 3.0*z1 + 2.0*z2)

    println("The optimization problem to be solved is:")
    println(model) # Shows the model constructed in a human-readable form

    status = optimize!(model)
    println("Objective value: ", objective_value(model)) 
    println("z1:",value(z1))
    println("z2:",value(z2))
end

main()

 

f:id:meison_amsl:20170801070633p:plain

function main3()
    println("JuMP sample5")
    model = Model(solver=MosekSolver())

    @variable(model, z1>=0.0)
    @variable(model, z2>=0.0)
    @objective(model, Min, z1-z2)
    @constraint(model, z1-z2>=30)
    @constraint(model, 2.0*z1+z2>=1)

    println("The optimization problem to be solved is:")
    println(model) # Shows the model constructed in a human-readable form

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(z2))
end

main3()

 

f:id:meison_amsl:20170801071107p:plain

function main5()
    println("JuMP sample7")
    model = Model(solver=MosekSolver())

    @variable(model, z1>=1.0)
    @variable(model, z2>=1.0)
    @objective(model, Min, z1^2+z2^2)

    println("The optimization problem to be solved is:")
    println(model) # Shows the model constructed in a human-readable form

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(z2))
end

 

f:id:meison_amsl:20170801072757p:plain

function main8()
    println("JuMP sample8")
    model = Model(solver=MosekSolver())

    @variable(model, z1>=0.0)
    @variable(model, z2>=0.0)
    @variable(model, z3>=0.0)
    @objective(model, Min, 0.5*(z1^2+z2^2+0.1*z3^2)+0.55*z3)
    @constraint(model, z1+z2+z3==1)

    println("The optimization problem to be solved is:")
    println(model) # Shows the model constructed in a human-readable form

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(z2))
    println("z3:",getvalue(z3))
end

 

JuMPによる最小二乗法解法サンプルコード

以前、cvxpyで解いた最小二乗法問題をJumpで解いてみました。

myenigma.hatenablog.com

 

using JuMP
using Mosek

function main()
    println("JuMP sample1")
    model = Model(solver=MosekSolver())
    m = 10
    n = 5
    A = rand(m, n)
    b = rand(m, 1)
    println(A)
    println(b)

    @variable(model, 0.0<=x[1:size(A,2)]<=1.0)
    @objective(model, Min, sum((A*x-b).^2))

    println("The optimization problem to be solved is:")
    println(model) # Shows the model constructed in a human-readable form

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("x:",getvalue(x))
end

main()

github.com

すると下記のように最小二乗法の計算結果が得られます。

modelをprintすると、

最適化問題が数式みたいに表示されるのはいい感じですね。

ただ行列をprintした時に、

行と列で綺麗に表示されないのが悲しい。。。

f:id:meison_amsl:20170730084633p:plain

f:id:meison_amsl:20170730084746p:plain

 

JuMPによる混合整数最適化問題の解法

続いて、

下記の2つの混合整数(Mixed integer)最適化問題をJuMPを使って解いてみます。

f:id:meison_amsl:20170803100420p:plain

f:id:meison_amsl:20170803100435p:plain

 

下記がJuliaのサンプルコードです。

今回はMixed integerを解くことができる。

MosekとCPLEXをソルバーとして使っています。

using JuMP
using Mosek
using CPLEX

function main1()
    println("JuMP mix integer sample1")
    model = Model(solver=MosekSolver())

    @variable(model, z1, Int)
    @variable(model, z2, Int)
    @objective(model, Min, -6*z1-5*z2)
    @constraint(model, z1+4*z2<=16)
    @constraint(model, 6*z1+4*z2<=28)
    @constraint(model, 2*z1-5*z2<=6)
    @constraint(model, 0<=z1<=10)
    @constraint(model, 0<=z2<=10)

    println("The optimization problem to be solved is:")
    println(model) # Shows the model constructed in a human-readable form

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(z2))
end

function main2()
    println("JuMP mix integer sample2")
    model = Model(solver=CplexSolver())

    @variable(model, z1>=0)
    @variable(model, z2>=0)
    @variable(model, z3, Bin)
    @objective(model, Min, -z1-2*z2)
    @constraint(model, 3*z1+4*z2-12<=4*z3)
    @constraint(model, 4*z1+3*z2-12<=4*(1-z3))

    println("The optimization problem to be solved is:")
    println(model)

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(z2))
    println("z3:",getvalue(z3))
end

main1()
main2()

github.com

 

JuMPによるナップサック問題の解法サンプルコード

以前cvxpyで解いたナップサック問題を

JuMPで解いてみました。

myenigma.hatenablog.com

 

using JuMP
using Mosek

function main()
    println("JuMP knapsack")
    model = Model(solver=MosekSolver())

    size = [21, 11, 15, 9, 34, 25, 41, 52]
    weight = [22, 12, 16, 10, 35, 26, 42, 53]
    capacity = 100

    @variable(model, x[1:length(size)], Bin)
    @objective(model, Max, dot(weight,x))
    @constraint(model, dot(size, x) <= capacity)

    println("The optimization problem to be solved is:")
    println(model) # Shows the model constructed in a human-readable form

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("Solution is:")
    for i = 1:length(size)
        println("x[$i] = ", getvalue(x[i]))
    end
end

main()

github.com

 

下記のようにちゃんとcvxpyの時と同じ結果が得られています。

f:id:meison_amsl:20170731015249p:plain

 

注意点としては、最適化変数をboolenにしたい時は、

Boolを与えるのではなく、JuMP専用のbinを型として与える必要があります。

 

続いて、同じくcvxpyでも実施した

一つの商品を複数回選んで良い場合のコードはこちらです。

using JuMP
using CPLEX

function main()
    println("JuMP knapsack")
    model = Model(solver=CplexSolver())

    size = [21, 11, 15, 9, 34, 25, 41, 52]
    weight = [22, 12, 16, 10, 35, 26, 42, 53]
    capacity = 100

    @variable(model, x[1:length(size)], Int)
    @constraint(model, x .>= 0)
    @objective(model, Max, dot(weight,x))
    @constraint(model, dot(size, x) <= capacity)

    println("The optimization problem to be solved is:")
    println(model) # Shows the model constructed in a human-readable form

    status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("Solution is:")
    for i = 1:length(size)
        println("x[$i] = ", getvalue(x[i]))
    end
end

main()

github.com

今回はなんとなくソルバーとしてCPLEXを使ってみました。

ソルバーオブジェクトの名前がCplexSolverであることに注意が必要です。

以前の結果と同じ結果が得られています。

f:id:meison_amsl:20170731021343p:plain

 

ナップサック問題におけるcvxpyとJuMPの計算速度比較

Juliaは計算が早いことで有名な言語なので、

cvxpyとJuMPで計算速度のベンチマークを取ってみました。

下記がcvxpyのコード

#
# A sample code to solve a knapsack_problem with cvxpy
#
# Author: Atsushi Sakai (@Atsushi_twi)
#

import cvxpy
import numpy as np
import random
import time

N_item = 30000

size = np.array([random.randint(1, 50) for i in range(N_item)])
weight = np.array([random.randint(1, 50) for i in range(N_item)])
capacity = 1000

x = cvxpy.Bool(size.shape[0])
objective = cvxpy.Maximize(weight * x)
constraints = [capacity >= size * x]

start = time.time()
prob = cvxpy.Problem(objective, constraints)
prob.solve(solver=cvxpy.ECOS_BB)
elapsed_time = time.time() - start
print("elapsed_time:{0}".format(elapsed_time) + "[sec]")
result = [round(ix[0, 0]) for ix in x.value]

print("status:", prob.status)
print("optimal value", prob.value)

下記がjuliaのコードです。

using JuMP
using CPLEX

function main()
    println("JuMP knapsack")
    model = Model(solver=CplexSolver())

    nitem = 30000
    capacity = 1000

    size = rand(1:50, nitem)
    weight = rand(1:50, nitem)

    @variable(model, x[1:length(size)], Bin)
    @objective(model, Max, dot(weight,x))
    @constraint(model, dot(size, x) <= capacity)

    @time status = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
end

main()

github.com

 

ナップサック問題の商品の数を30000個にした場合の比較です。

cvxpy JuMP
計算時間 14.2[s] 1.4[s]

十倍以上早いですね。。。

Juliaのコンパイル時間をいれても3秒ほどだったので

かなり早くなっています。

CPLEXがすごいという噂もありますが。

 

最大流問題(Max Flow Problem)の解き方

有名な線形計画問題に最大流問題がありますが、

こちらの解き方に関しては、下記の記事を参照ください。

myenigma.hatenablog.com

 

最小コストネットワークフロー問題

有名な線形計画問題に、電力網の最適化などにに使われる

最小コストネットワークフロー問題がありますが、

こちらの解き方に関しては、下記の記事を参照ください。

myenigma.hatenablog.com

 

より詳しく線形計画法を学びたい人は

下記の書籍がおすすめです。

 

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

 

MyEnigma Supporters

もしこの記事が参考になり、

ブログをサポートしたいと思われた方は、

こちらからよろしくお願いします。

myenigma.hatenablog.com