PCガジェットはいつでも安いAmazonで!

Calcのソルバーをプログラムで繰り返す【初期値を変えて複数解を自動探索】

スタッフ厳選!日本HP おすすめパソコン

icon icon

今回はLibreOffice Calcの非線形ソルバーをプログラムで繰り返して複数解を自動探索してみるという内容です。

目次

マクロで自動的にソルバーを回したい

以前複数解を自動探索するプログラムをVBAとMicrosoft Excelの組み合わせで作ってみました。

あわせて読みたい
ソルバーをVBAで繰り返す【初期値を変えて複数解を自動探索】 今回はエクセルのソルバーをVBAで操作して、繰り返し処理でソルバーを繰り返し実行し、繰り返しごとに初期値を変えて複数解を自動探索するようにします。 最終的にはコ...

これをLibreOfficeのCalcのソルバーでできないかなというのが今回の内容です。

なお通常のダイアログから普通にソルバーを実行したい場合は以下の記事をご覧ください。

あわせて読みたい
Calcのソルバーでどこまでできるか実験 はじめに 今回はまたソルバーネタ。 LibreOfficeのCalcにもソルバーがあるので、これがどのくらいの性能なのか検証してみようという回。 Calcは最近のアップデートでDEP...

まずは簡単な問題でプログラムからソルバーを動かしてみる

Calcのソルバーをプログラムから動かすという話をわかりやすく体験するために、簡単な問題をやってみます。

\[f=x^2\]

\[-10\le x \le 10\]

答えは\(x_G =0\)ですね。

これをCalcの非線形ソルバーであるDEPSでプログラムで求めてみます。

プログラム

まずシート(シート名「Sheet1」)に以下のように打ち込みます。

B3:6

C3:=B3*B3

そして以下のプログラムを「ツール」→「マクロ」→「マクロの編集」をクリックして出てくるマクロ入力画面に入力します。

Option Explicit

Sub SolverMacro
' Part1
  Dim orSolver As Object
  Dim orDocument As Object
  Dim orSheet As Object
  Dim i As Integer

  Dim orObjectiveCell As Object

' Part2
  Dim orVariableCells(0) as new com.sun.star.table.CellAddress
  Dim orConstraints(1) as new com.sun.star.sheet.SolverConstraint

  Set orSolver = CreateUnoService("com.sun.star.comp.Calc.NLPSolver.DEPSSolverImpl")

  orDocument = ThisComponent
  orSolver.Document = orDocument

  orSheet = orDocument.Sheets.getByName("Sheet1")

' Part3
  orObjectiveCell =  orSheet.getCellByPosition(2,2).CellAddress
  orSolver.Objective = orObjectiveCell

  orSolver.Maximize = false
  orVariableCells(0) = orSheet.getCellByPosition(1,2).CellAddress
  orSolver.Variables = orVariableCells

  orConstraints(0).Left = orSheet.getCellByPosition(1,2).CellAddress
  orConstraints(0).Operator = com.sun.star.sheet.SolverConstraintOperator.LESS_EQUAL
  orConstraints(0).Right = 10

  orConstraints(1).Left = orSheet.getCellByPosition(1,2).CellAddress
  orConstraints(1).Operator = com.sun.star.sheet.SolverConstraintOperator.GREATER_EQUAL
  orConstraints(1).Right = -10

  orSolver.Constraints = orConstraints()

' Part4
  orSolver.AgentSwitchRate = 1.0
  orSolver.DECR = 0.9

  orSolver.EnhancedSolverStatus=false
  orSolver.solve

' Part5
  If orSolver.Success then
  	for i = 0 to UBound(orSolver.Variables)
  		orSheet.getCellByPosition(orSolver.Variables(i).Column,orSolver.Variables(i).Row).Value = orSolver.Solution(i)
  	next
  Else
     MsgBox "No solution was found."
  End If

End Sub

Part1

  Dim orSolver As Object
  Dim orDocument As Object
  Dim orSheet As Object
  Dim i As Integer

  Dim orObjectiveCell As Object

各種オブジェクトの入る変数を定義します。

orSolverがソルバーを計算する本体のオブジェクト、orDocumentとorSheetがファイルやシート情報などを格納するオブジェクトです。

orObjectiveCellに目的関数の値が入るセルを格納します。

Part2

  Dim orVariableCells(0) as new com.sun.star.table.CellAddress
  Dim orConstraints(1) as new com.sun.star.sheet.SolverConstraint

  Set orSolver = CreateUnoService("com.sun.star.comp.Calc.NLPSolver.DEPSSolverImpl")

  orDocument = ThisComponent
  orSolver.Document = orDocument

  orSheet = orDocument.Sheets.getByName("Sheet1")

今回は変数が一つなのでoVariableCells(0)となっていますが、これが3変数ならoVariableCells(2)、5変数ならoVariableCells(4)としてください。

orConstraints(1)も同様で、制約条件が2つなら1を、5つなら4を入れてください。

今回は以下の条件をソルバーに設定する必要があるので、制約条件は2つです。

\[-10 \le x\]

\[10 \ge x\]

com.sun.star.table.CellAddressは変数につけるおまじない、com.sun.star.sheet.SolverConstraintは制約条件につけるおまじないと考えましょう。非線形ソルバーのプログラムから必要なクラスとかを呼び出しているだけです。

com.sun.star.comp.Calc.NLPSolver.DEPSSolverImplはDEPSを使う時のおまじないです。

他に線型計画法ソルバーなら「com.sun.star.comp.Calc.LpsolveSolver」、SCOソルバーなら「com.sun.star.comp.Calc.NLPSolver.SCOSolverImpl」を入力します。

「orDocument = ThisComponent」の部分は現在のファイルを参照するという意味です。

「orDocument.Sheets.getByName(“Sheet1”)」の部分は使用するシートの名前を記述します。今回は「Sheet1」で実行するので「Sheet1」と記述します。

Part3

  orObjectiveCell =  orSheet.getCellByPosition(2,2).CellAddress
  orSolver.Objective = orObjectiveCell

  orSolver.Maximize = false
  orVariableCells(0) = orSheet.getCellByPosition(1,2).CellAddress
  orSolver.Variables = orVariableCells

  orConstraints(0).Left = orSheet.getCellByPosition(1,2).CellAddress
  orConstraints(0).Operator = com.sun.star.sheet.SolverConstraintOperator.LESS_EQUAL
  orConstraints(0).Right = 10

  orConstraints(1).Left = orSheet.getCellByPosition(1,2).CellAddress
  orConstraints(1).Operator = com.sun.star.sheet.SolverConstraintOperator.GREATER_EQUAL
  orConstraints(1).Right = -10

  orSolver.Constraints = orConstraints()

「orObjectiveCell = orSheet.getCellByPosition(2,2).CellAddress」は目的関数の値のセルの情報を入力します。今回はC3セルなので(2,2)となっています。

「orSolver.Maximize = false」は目的関数を最小化するという意味です。trueにすると最大化になります。

「orVariableCells(0) = orSheet.getCellByPosition(1,2).CellAddress」は変数セルを入力します。今回はB3セルなので(1,2)となっています。列が最初の数値、行が後の数値になっています。

もし変数セルが複数あるなら「orVariableCells(1) = orSheet.getCellByPosition(1,3).CellAddress」(B4セル)などをいちいち指定する必要があります。

「orConstraints(0).Left = orSheet.getCellByPosition(1,2).CellAddress」に関してはソルバーの手動設定ダイアログの制約条件の入力を考えてみるとわかりやすいでしょう。

左側にB3セルと指定します。

そして不等号などの入力を「orConstraints(0).Operator = com.sun.star.sheet.SolverConstraintOperator.LESS_EQUAL」としています。

「小なり」が「com.sun.star.sheet.SolverConstraintOperator.LESS_EQUAL」、「等号」が「com.sun.star.sheet.SolverConstraintOperator.EQUAL」、「大なり」が「com.sun.star.sheet.SolverConstraintOperator.GREATER_EQUAL」、「整数」が「com.sun.star.sheet.SolverConstraintOperator.INTEGER」、「バイナリ」が「com.sun.star.sheet.SolverConstraintOperator.BINARY」となります。

「orConstraints(0).Right = 10」で制約条件の右側を設定します。今回は10以下を表すので10が入力されています。

2つ目の制約条件も同様に入力します。

Part4

  orSolver.AgentSwitchRate = 1.0
  orSolver.DECR = 0.9

  orSolver.EnhancedSolverStatus=false
  orSolver.solve

DEPSのオプションの設定です。詳細は以下のページに全ての設定項目があるのでご確認ください。

参考(1):Apache OpenOffice Wiki, NLPSolver

「orSolver.AgentSwitchRate」が「エージェントスイッチレート」つまり「DEの確率」。100%をDEで使いたいので1.0を入力しています。

「orSolver.DECR = 0.9」交叉率CRの値です。その他に「orSolver.DEFactor = 0.5」なんて設定もできます。これはスケーリング係数というものです。

交叉率とスケーリング係数(スケール係数)の意味とうまくいきやすいパラメータ値に関しては以下の文献をご確認ください。今回は文献を参考にしてスケーリング係数はデフォルト、交叉率は0.9としています。ここをいじりたいときに参考にしてみてください。

参考(2):田川 聖治, 距離に基づく生存選択を用いたDifferential Evolutionの構成法, 電気学会論文誌C(電子・情報・システム部門誌)/130 巻 (2010) 5 号

また意外に重要なのが「orSolver.EnhancedSolverStatus=false」です。これがデフォルトだとtrueになっていて、ソルバーが解決するとダイアログが出てきてしまいます。連続でソルバーを回す時は出てきてほしくないので、これをfalseとすることでダイアログが表示されないようになります。

「orSolver.solve」でソルバーを実行します。

Part5

  If orSolver.Success then
  	for i = 0 to UBound(orSolver.Variables)
  		orSheet.getCellByPosition(orSolver.Variables(i).Column,orSolver.Variables(i).Row).Value = orSolver.Solution(i)
  	next
  Else
     MsgBox "No solution was found."
  End If

もし解が求まったなら「orSolver.Success」はtrueになっているのでこのIf文が実行されます。

「UBound(orSolver.Variables)」が変数の数を表します。

「orSheet.getCellByPosition(orSolver.Variables(i).Column,orSolver.Variables(i).Row).Value = orSolver.Solution(i)」この左辺がようするに変数セルです。右辺がソルバーが解決して得られた解を表します。変数セルが縦に並んでいるか横に並んでいるかで書き方がちょっと変わるので適宜修正してください。

oSolver.Objectiveには自動的に求めた解の目的関数の値が入るようです。

参考サイト一覧

今回の記事で参考にしたサイトを以下にまとめておきます。

まあ色々なパラメーターはソルバーオブジェクト(例えば今回は変数orSolver)に「.〜」と書けば動くみたいです。(〜に各種パラメーターが入る)

ソルバーをプログラムで繰り返してみる

あわせて読みたい
ソルバーをVBAで繰り返す【初期値を変えて複数解を自動探索】 今回はエクセルのソルバーをVBAで操作して、繰り返し処理でソルバーを繰り返し実行し、繰り返しごとに初期値を変えて複数解を自動探索するようにします。 最終的にはコ...

上の記事でやったことをCalcのソルバーでやってみます。

  • B3:0
  • C3:=COS(B3*PI()/180)

コサインカーブで複数解を求める問題です。

  • 目的関数:\(y=\cos x\)→最小化
  • 制約条件:\(0 \leq x \leq 1080\)
  • 変数:\(x\)

これを0≦x≦1080で求めるプログラムが以下となります。

Option Explicit

Sub SolverMacro
  Dim orSolver As Object
  Dim orDocument As Object
  Dim orSheet As Object
  Dim i As Integer
 
  Dim x(100) As Double
  Dim f(100) As Double

  Dim orObjectiveCell As Object

  Dim orVariableCells(0) as new com.sun.star.table.CellAddress
  Dim orConstraints(1) as new com.sun.star.sheet.SolverConstraint

  Set orSolver = CreateUnoService("com.sun.star.comp.Calc.NLPSolver.DEPSSolverImpl")

  orDocument = ThisComponent
  orSolver.Document = orDocument

  orSheet = orDocument.Sheets.getByName("Sheet1")

  orObjectiveCell =  orSheet.getCellByPosition(2,2).CellAddress
  orSolver.Objective = orObjectiveCell

  orSolver.Maximize = false
  orVariableCells(0) = orSheet.getCellByPosition(1,2).CellAddress
  orSolver.Variables = orVariableCells
  
  orSolver.AgentSwitchRate = 1.0
  orSolver.DECR = 0.9
  
  for i = 0 to 12
    orConstraints(0).Left = orSheet.getCellByPosition(1,2).CellAddress
    orConstraints(0).Operator = com.sun.star.sheet.SolverConstraintOperator.LESS_EQUAL
    orConstraints(0).Right = (i + 1) * 90

    orConstraints(1).Left = orSheet.getCellByPosition(1,2).CellAddress
    orConstraints(1).Operator = com.sun.star.sheet.SolverConstraintOperator.GREATER_EQUAL
    orConstraints(1).Right = i * 90

    orSolver.Constraints = orConstraints()

    orSolver.EnhancedSolverStatus=false
    orSolver.solve
    x(i) = orSolver.Solution(0)
    f(i) = orSolver.ResultValue
  next
  If orSolver.Success then
  	for i = 0 to 12
  		orSheet.getCellByPosition(1, i + 4).Value = x(i)
  		orSheet.getCellByPosition(2, i + 4).Value = f(i)
  	next
  Else
     MsgBox "No solution was found."
  End If

End Sub

上のCalcの基本的な解説プログラムとの違いだけ解説します。

  for i = 0 to 12
    orConstraints(0).Left = orSheet.getCellByPosition(1,2).CellAddress
    orConstraints(0).Operator = com.sun.star.sheet.SolverConstraintOperator.LESS_EQUAL
    orConstraints(0).Right = (i + 1) * 90

    orConstraints(1).Left = orSheet.getCellByPosition(1,2).CellAddress
    orConstraints(1).Operator = com.sun.star.sheet.SolverConstraintOperator.GREATER_EQUAL
    orConstraints(1).Right = i * 90

    orSolver.Constraints = orConstraints()

    orSolver.EnhancedSolverStatus=false
    orSolver.solve
    x(i) = orSolver.Solution(0)
    f(i) = orSolver.ResultValue
  next

「orConstraints(0).Right = (i + 1) * 90」と「orConstraints(1).Right = i * 90」が変数の右側と左側の限界値を表します。90°刻みでソルバーを回すということです。

x(i)とf(i)に求めた解を格納していきます。今回は1変数だけなのでorSolver.Solution(0)でいいですが、複数の変数がある場合はorSolver.Solution(1)などと別途格納する必要があります。

orSolver.ResultValueは解決で求まった解のときの目的関数の値が入っているのでそれをf(i)に格納します。

  If orSolver.Success then
  	for i = 0 to 12
  		orSheet.getCellByPosition(1, i + 4).Value = x(i)
  		orSheet.getCellByPosition(2, i + 4).Value = f(i)
  	next
  Else
     MsgBox "No solution was found."
  End If

上の部分はあまり説明する必要は無いと思いますが、やっていることとしてはB列とC列の下の方に解と目的関数の値を出力しているだけです。

結果は以下となります。

ちゃんと180°と540°と900°で最適解のときの目的関数値である-1が求められていますね。

後は適宜ご自身の最適化問題で調整してみましょう。

まとめ【Calcのソルバーでもまあまあ使える】

今回はLibreOffice Calcの非線形ソルバーをプログラムで動かしてみるという内容でした。

もし最適解を何回も自動的に求めたいというときに参考になりましたら幸いです。

目次