【Julia言語入門】Juliaとmodelicaで魚の体重をPI制御する

今回もこれまでに引き続いてJuliaとmodelicaの連携を実践してみます。

今回は「魚の体重」をJuliaとmodelicaを使ってPI制御してみます。

魚の体重の数理モデル

魚の体重の増加量を表す数理モデルとしてフォン・ベルタンフィーモデルが知られています。

$\frac{dw}{dt} = \alpha w^{2/3} - \beta w$

ここで $w$ は魚の体重を示します。また右辺第1項は栄養分による体重の増加を、第2項は呼吸による体重ロスを表します。

今回はエサの量≒係数 $\alpha$ を入力とし、魚の体重 $w$ をコントロールします。

過去に同様の数理モデルを使ってフィードバック制御をした際の記事は以下になります。この記事ではJuliaのみでのシミュレーションでした。

otepipi.hatenablog.com

システムの全体像

今回構築するシステムの全体像は以下のようになります。

f:id:Otepipi:20210209064303p:plain

modelicaで以下の数理モデルを構築します。オブジェクトとして『魚』『PIコントローラー』『給餌器』の３つがあります。それぞれの役割は以下です。

魚……餌量と体重の関係を実装します。フォン・ベルタンフィーモデルを実装します。
PIコントローラー……魚の体重が目標値になるように、餌量を給餌器に指示します。
給餌器……PIコントローラーからの指示にあった餌量を、魚に供給します。

そして、上記modelicaの数理モデルを、Juliaから操作します。魚の体重の目標値を与えたり、実行司令や計算結果の受け渡し、計算結果の描画等を実行します。

modelicaのモデル

今回作成したmodelicaのモデルを紹介します。なお、コードは以下書籍を参考にしています。

Modelicaによるシステムシミュレーション入門 -モデルベース開発のための物理システムモデリング- (MBD Lab Series)

作者:Peter Fritzson
発売日: 2015/07/15
メディア: 単行本（ソフトカバー）

魚

フォン・ベルタンフィーモデルを実装しています。魚の体重の初期値は10としています。

model Fish
        WeightSignal wSensor ;
        Foodflow foodIn;
        parameter Real beta = 1;
        Real w(start=10)"Fish weight";
    equation
        der(w) = foodIn.fflow*w^(2/3) - beta*w;
        wSensor.val = w;
    end Fish;

PIコントローラー

PI制御を実装しています。今回比例ゲイン=0.01, 積分時間=1 としています。

model PIcountinuousController
        WeightSignal cIn;
        FoodSignal cOut;
        parameter Real K=0.01 "Gain";
        parameter Real T=1 "Time Constans";
        Real ref "Reference level";
        Real error "Deviation level";
        Real outCtr "Output contril signal";
        Real x;//State variable of continuous PI controller
    equation
        error = ref - cIn.val;
        der(x) = error/T;
        outCtr = K*(error +x);
        cOut.fval = outCtr;
    end PIcountinuousController;

給餌器

PIコントローラーからの餌量司令によって餌を供給します。最小の餌量を0, 最大の餌量を100とするための関数LimitValueを実装しています。

model FoodSource
        Foodflow foodOut;
        FoodSignal order;
        parameter Real foodmin = 0;
        parameter Real foodmax = 100;
    equation
        foodOut.fflow = LimitValue(foodmin,foodmax,order.fval);
    end FoodSource;

    function LimitValue
        input Real pMin;
        input Real pMax;
        input Real p;
        output Real pLim;
    algorithm
        pLim := if p>pMax then pMax
                else if p<pMin then pMin
                else p;    
    end LimitValue;

コネクタ類

今回、『魚』『PIコントローラー』『給餌器』のオブジェクト間を、『魚の体重情報』『餌』『餌量の司令値情報』の3つが行き交いますのでそれぞれのコネクタを実装してやります。

connector WeightSignal"Reading fish weight"
        Real val;
    end WeightSignal;

    connector FoodSignal"Signal to food mass"
        Real fval;
    end FoodSignal;

    connector Foodflow "Food flow"
        Real fflow;
    end Foodflow;

まとめのモデル

最後に各オブジェクトを配置したモデルFishPIを作ります。このモデルを実行することでシミュレーションが可能です。体重の目標値m_iをinputとしてJuliaから入力するようにします。

 model FishPI
        Fish fish;
        FoodSource source;
        PIcountinuousController pi;
        input Real m_i;
        output Real wei;
    equation
        connect(source.foodOut, fish.foodIn);
        connect(fish.wSensor, pi.cIn);
        connect(pi.cOut, source.order);
        wei = fish.w;
        pi.ref = m_i;
    end FishPI;

Juliaの実装とシミュレーションの実行

今回もopen modelicaのJulia APIであるOMJuliaを使用して、Juliaからmodelicaモデルを実行します。これまでの関連記事は以下です。

otepipi.hatenablog.com

コードは以下のようになります

using Plots
using OMJulia


tnk = OMJulia.OMCSession();
tnk.ModelicaSystem("FishWeight.mo","FishWeight.FishPI");

tnk.setSimulationOptions(["stopTime=100"]);
tnk.setInputs("m_i = [(0,80),(20,80),(50,40),(80,40)]")
tnk.simulate();
tm,wei,ref = tnk.getSolutions(["time","wei","m_i"]);

plot(tm,[wei,ref],title="PI control",xlabel="time[s]",
    label = ["weight" "set point"],lw=3)