読者です 読者をやめる 読者になる 読者になる

MyEnigma

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

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

optimization Programming Python

文科系の線形計画法入門

文科系の線形計画法入門

目次

はじめに

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

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

myenigma.hatenablog.com

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

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

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

 

線形計画法の概要

線形計画法は、

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

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

f:id:meison_amsl:20161025051050p:plain

一つ注意点としては、変数xiは0以上であるという制約があるということです。

 

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

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

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

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

f:id:meison_amsl:20161203063324p:plain

 

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

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

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

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

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

 

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

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

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

最適解がある

ということです。

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

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

 

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

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

シンプレックス法などが有名です。

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

 

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

 

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

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

これなら分かる最適化数学―基礎原理から計算手法まで

これなら分かる最適化数学―基礎原理から計算手法まで

文科系の線形計画法入門

文科系の線形計画法入門

 

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com