1. はじめに

IvyFEMで電磁導波路の伝送特性を計算できるようにしたので、せっかく有限要素法を使っているんだから最適化に使ってみたいと思い試してみました。

2. 最適化問題

直角コーナーベンド導波路の反射電力が対象周波数領域で最も低くなるような曲がり部の形状を求める。
曲がり部を等角度で図のように分割することを考える。このとき、内側コーナー点から外側コーナーの境界までの長さd1, d2, ... , dnを定義できる。(図ではn = 4)
この{d}の集合を変数とし、
　　v({d}) = ∑f=f1,fp {|S11|^2 }　　　　　(1)
　　f1, ..., fp : 周波数ポイント
|S11| ポート1の反射係数の絶対値
|S21| ポート2 → 1の透過係数の絶対値
|S11|^2 + |S21|^2 = 1
が0(最小)になるように{d}を求めればそれが最適形状である。

2. Microsoft.SolverFoundationのインストール

NuGetパッケージとしてインストールできます。

3. Microsoft.SolverFoundationを使った最適化関数の実装

NonlinerOptimizingProblem.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using Microsoft.SolverFoundation.Solvers;
using Microsoft.SolverFoundation.Services;

namespace IvyFEMProtoApp
{
    class NonlinerOptimizingProblem
    {
        public static void Solve(
            int paramCnt, double[] initialValues, double[] minValues, double[] maxValues, int maxIter,
            Func<INonlinearModel, int, ValuesByIndex, bool, object, double> callbackFunc, object callbackArg)
        {
            System.Diagnostics.Debug.Assert(paramCnt == initialValues.Length);

            NelderMeadSolver solver = new NelderMeadSolver();

            // 目的変数
            int objVId;
            solver.AddRow("obj", out objVId);
            solver.AddGoal(objVId, 0, true); // true: minimize

            for (int i = 0; i < paramCnt; i++)
            {
                // 変数
                int vId;
                solver.AddVariable("v" + i, out vId);

                // 範囲
                solver.SetBounds(vId, minValues[i], maxValues[i]);

                // 初期値
                solver.SetValue(vId, initialValues[i]);
            }

            // 目的関数をセット
            Func<INonlinearModel, int, ValuesByIndex, bool, double> callbackFunc2 =
                (model, rowVId, values, newValues) =>
                {
                    return callbackFunc(model, rowVId, values, newValues, callbackArg);
                };
            solver.FunctionEvaluator = callbackFunc2;

            // 中止関数
            Func<bool> abortFunc = () =>
            {
                return solver.IterationCount >= maxIter;
            };

            // 解く
            var solverParams = new NelderMeadSolverParams(abortFunc);
            INonlinearSolution solution = solver.Solve(solverParams);

            {
                System.Diagnostics.Debug.WriteLine("The Result is " + solution.Result + ".");
                System.Diagnostics.Debug.WriteLine("The minimium objective value is " +
                                  solution.GetValue(objVId) + ".");
                for (int i = 0; i < paramCnt; i++)
                {
                    string vKey = "v" + i;
                    int vId = solver.GetIndexFromKey(vKey);
                    System.Diagnostics.Debug.WriteLine(vKey + " = " +
                        solution.GetValue(vId));
                }
            }
        }
    }
}

4. 最適化関数の呼び出し

Problem.OptimizeProblem1.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using IvyFEM;
using IvyFEMProtoApp;
using OxyPlot;
using OxyPlot.Axes;
using OxyPlot.Series;
using Microsoft.SolverFoundation.Solvers;
using Microsoft.SolverFoundation.Services;

namespace IvyFEMProtoApp
{
    partial class Problem
    {
        public void Optimize1Problem(MainWindow mainWindow)
        {
            int paramCnt = 5;
            double[] ds = new double[paramCnt];
            double[] minds = new double[paramCnt];
            double[] maxds = new double[paramCnt]; 
            for (int i = 0; i < paramCnt; i++)
            {
                // 直角コーナー
                double theta = (i + 1) / (double)paramCnt * Math.PI / 4.0;
                double x = 1.0;
                double y = x * Math.Tan(theta);
                ds[i] = Math.Sqrt(x * x + y * y);
                minds[i] = 0.1;
                maxds[i] = 2.0;
            }
            double totalS11;
            WaveguideProblemToBeOptimized(mainWindow, paramCnt, ds, out totalS11);

            int maxItr = 50;
            NonlinerOptimizingProblem.Solve(paramCnt, ds, minds, maxds, maxItr, ValueFunctionToBeOptimized, mainWindow);
        }

        private double ValueFunctionToBeOptimized(INonlinearModel model, int rowVId,
                ValuesByIndex values, bool newValues, object arg)
        {
            MainWindow mainWindow = arg as MainWindow;

            System.Diagnostics.Debug.Assert(model.GetKeyFromIndex(rowVId) as string == "obj");
            int paramCnt = model.VariableCount;
            double[] ds = new double[paramCnt];
            for (int i = 0; i < paramCnt; i++)
            {
                string vKey = "v" + i;
                int vId = model.GetIndexFromKey(vKey);
                ds[i] = values[vId];
            }

            double totalS11;
            WaveguideProblemToBeOptimized(mainWindow, paramCnt, ds, out totalS11);
            return totalS11;
        }