MyEnigma

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

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

文科系の線形計画法入門

文科系の線形計画法入門

目次

はじめに

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

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

myenigma.hatenablog.com

今回はもっと多くの分野で使われている最適化問題の一つである

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

同じくPythonでの解き方を紹介したいと思います。

 

線形計画法の概要

線形計画法は、

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

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

f:id:meison_amsl:20161025051050p:plain

 

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

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

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

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

f:id:meison_amsl:20161203063324p:plain

 

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

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

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

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

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

 

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

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

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

最適解がある

ということです。

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

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

 

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

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

シンプレックス法や

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

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

myenigma.hatenablog.com

 

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

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

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

 

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

先日、二次計画法の際に紹介した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が等号制約を表しています。

 

#! /usr/bin/python 
# -*- coding: utf-8 -*- 
u""" 
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(sol["x"])
print(sol["primal objective"])

 

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

f:id:meison_amsl:20161025092259p:plain

先程の書籍にある通り、

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

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

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

 

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

github.com

 

Python製線形計画法モデリングライブラリPuLP

上記のcvxoptを使えば、線形計画法を解くことができますが、

解きたい問題をわざわざ上記の線形計画法の標準形に

変形する必要があります。

 

そんな時に便利なのが、

線形計画法PythonモデリングライブラリPuLPです。

このPuLPを使うことで、

制約条件の式や目的関数の数式をそのまま利用して、

線形計画法を解くことができます。

 

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 Mosek

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

    @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 = solve(model)
    println("Objective value: ", getobjectivevalue(model)) 
    println("z1:",getvalue(z1))
    println("z2:",getvalue(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 Supporters

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

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

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

myenigma.hatenablog.com