複数の方程式からなる​関数を用いたフィッテ​ィング方法

12 views (last 30 days)
Hiroki ENDO
Hiroki ENDO on 13 Jul 2020
Commented: Hiroki ENDO on 14 Jul 2020
ある実測値に対して,複数の方程式(連立方程式)を用いてフィッティングを行うようにプログラミングしています.しかし,エラーが出て,なかなか解決できていません.エラーは「入力因数が不足している」とのことです.以下のプログラムを正常にするには,どのようにすれば良いでしょうか?
 ・関数1(目的関数):予備関数1(微分方程式)からAを算出し,Aから予備関数2(一次関数)を用いてsolを出力
 ・関数2:関数1における一部のパラメータを最適化変数にして,関数1を最適化関数にする.
その後,Bと実測値Bから最小二乗法で残差を計算し,最適化問題を作成する.
そして,Fitting変数を求める.
 ・エラー:最適化関数に変化するときに,入力因数不足となる;入力因数を増やしたり,減らしたり,ほかの変数にしたりしたが変化なし
[関数1]
function [T, sol] = fun(tspan, re, Azero, L, c, lambda, Vol, Vism);
re = optimvar('re',4,'LowerBound',0,'UpperBound',inf);
C = re(1);
S = re(2);
M1 = re(3);
M2 = re(4);
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-9);
[T, Av] = ode45('Adot', tspan, tens2vec(Azero), options, L, C, c, lambda, S); %予備関数1 微分方程式
for i = 1:length(T);
[Tau(i,:), Eta(i,:)] = Lipscomb(Av(i,:), L, c,S, Vol, Vism, M1, M2); %予備関数2 線形方程式
end
sol = Eta(:,6);
end
[関数2]
function [C, S, a, Vsteady, M1, M2, Rt] = EFitting(tspan, Azero, L, c, lambda, Vol, Vism, Eta_exp);
type fun
re = optimvar('re',4,'LowerBound',0,'UpperBound',inf); フィッティングパラメータ
re(1).LowerBound = 0;
re(1).UpperBound = 1;
[T, sol] = fcn2optimexpr('fun', tspan, re, Azero, L, c, lambda, Vol, Vism); %関数1を最適化関数にする
for i = 1:length(T);
Eta_temp(i:1) = Eta_exp(i:3); %実測値の3列目を抽出
end
obj = sum(sum((sol - Eta_temp).^2)); 
prob = optimproblem("Objective",obj); 
re0.re = [C S]; 
[resol,sumsq] = solve(prob,re0); %最適化関数を解く
C= resol(1); 出力1
S= resol(2); 出力2
end

Answers (1)

Shojiro SHIBAYAMA
Shojiro SHIBAYAMA on 14 Jul 2020
この関数のドキュメント読むと,以下の利用例がありました.
myfcn = fcn2optimexpr(@RtoODE,r,tspan,y0);
これを踏まえると,
[T, sol] = fcn2optimexpr(@fun, tspan, re, Azero, L, c, lambda, Vol, Vism); %関数1を最適化関数にする
と書くと良いと思うのですが,いかがでしょうか.
また,MATLABエディタの左枠に行数横の短い横棒線をクリックすると,そこでプログラム実行を止めてデバッグが可能です.こちらも利用して問題発生箇所を特定すると良いと思います.
  1 Comment
Hiroki ENDO
Hiroki ENDO on 14 Jul 2020
ありがとうございます.
ご指摘の通りにいたしましたが,効果はありませんでした.
ODEの場合は,
[T,Av] = ode45('Adot', tspan, tens2vec(Azero), options, L, C, c, lambda, S);
で動作することは確認済みです.

Sign in to comment.

Categories

Find more on Optimization Toolbox in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!