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;
        }

WaveguideProblemToBeOptimized()は{d}に対して、(1)式を求めるメソッドです。IvyFEMを使って実装しました。(IvyFEMのサンプルにある直角コーナーベンドとほぼ同じコードです。)

5. 計算結果

曲がり部の形状と|S11|、|S21|周波数特性の計算結果を示しておきます。

5.1. 初期形状

5.1. 最適化反復途中1

5.2. 最適化反復途中2

5.3. 最適化反復途中3

5.4. 最適化反復途中4

5.5. 最適化反復途中5

5.6. 最適化反復最大(50回)

円弧になるかなと想像しましたが3つ共振器をつなげたような形になりました。

6. まとめ

Microsoft.SolverFoundationを使って直角コーナーベンド導波路の伝送特性の最適化形状を求めるプログラムを実装しました。

ソースコードを
こちら(github: IvyFEMProtoApp)
にあげました。