スポンサーリンク
はじめに
前回の記事の文末に、「”バネ定数を設計変数としたときの構造最適化の手法”について説明しようと思う」と書いてしまいましたが、先にロバスト最適化について説明しようと思います。（予定を変更して、ごめんなさい。）
なので、今回は「ロバスト最適化」について説明します。ロバスト最適化の手法はいろいろありますが、下記記事で説明した「MATLABで学ぶ振動工学 質量付加を利用した構造最適化」を参考に説明します。
前提となる理論の復習
（復習）質量付加による伝達関数の変化
構造体の点iに質量mが付加されたとき、力Fが作用している点hと点j間の伝達関数newGjhは質量付加前の伝達関数Gjhで表すことができます。
$$
newG_{jh}=G_{jh}-\frac{G_{ji}G_{ih}}{G_{ii}-\frac{1}{ω^2 m}}
$$
図1 質量付加による伝達関数の変化
（復習）構造最適化の概要
最適化問題で値を変更しても良い変数を「設計変数」と呼びます。下図の付加質量mを変更してよい場合、mが設計変数となります。
図2 例題
質量付加する前の伝達関数G_jhの付加質量mに対する感度を求めます。ただし、直接的に求められないので、冒頭で示した質量付加後の伝達関数newG_jhの式から求めます。
質量付加後の伝達関数newG_jhの付加質量mに対する感度は、伝達関数newG_jhを付加質量mの偏微分で表されます。従って以下になります。
$$
\frac{\partial newG_{jh}}{\partial m} = \frac{G_{ji}G_{ih}}{ω^2(m G_{ii}-\frac{1}{ω^2} )^2}
$$
このmが微小のとき（mが0に近いとき）が”伝達関数G_jhの付加質量mに対する感度”なので、以下で表現できます。
$$
\frac{\partial G_{jh}}{\partial m} = ω^2 G_{ji} G_{ih}
$$
さらに、mが微小のときはG_{jh}をテイラー展開して、1次までを考慮した近似式で予測できます。従って下式が成立します。
$$
newG_{jh} =G_{jh} + \frac{\partial G_{jh}}{\partial m} dm
$$
（初級編）ロバスト最適化の概要
ロバスト最適化の概要
本記事を読んでいる学生さんの中には下記のような状態の人もいるかと思います。
学生
そもそもロバストって何ですか？
この質問への回答は難しいですね。というのも、分野や論文によって定義が異なるからです。ただ、本サイトでは「誤差因子に対する”何か”の安定性」という意味です。”何か”とは振動振幅であったり、共振周波数であったりします。
学生
誤差因子って何ですか？
誤差因子は誤差を持つ要因だと思ってください。
例えば、機械構造物（車や飛行機など）は、機種が同じならば同じ性能のような気がしますが、1つ1つの機械構造物は各部品の製造誤差や組付け誤差によって性能が微妙に違います。この「性能のばらつき≒ロバスト」です。ロバスト性が高いというのは、機械的な特性や性能のばらつきが小さいという意味になります。そして、誤差因子というのは、製造誤差や組付け誤差を指します。
ただし、ロバスト最適化で扱えるのは、モデル化できる問題のみです。例えば、バネマスモデルであれば、mかkに反映させられる問題のみ扱えることになります。そのため、部品の製造誤差はmやkに簡単に反映することができますが、組付け誤差は何かしらの方法でmやkに落とし込む工夫が必要です。（具体的には、組付け時のねじ締結が緩い場合はkが小さいなど、組付け誤差とモデルの関係性を明確にしないといけません。）この工夫が難しいので、一般的には部品などの製造誤差（製造ばらつき）のみを扱うことが多いです。
ロバスト最適化の理論
繰り返しになりますが、ロバスト最適化の理論は多数あります。というより、前提とする物理式によって、展開する理論式が変わります。また、本記事では難しい内容はできるだけ省いてわかりやすく説明したいので、専門家にはツッコミどころがあるかもしれませんが、ご了承ください。
今回は下式（構造最適化に用いた式）をもとにロバスト最適化を考えます。
$$
newG_{jh} =G_{jh} + \frac{\partial G_{jh}}{\partial m} dm
$$
そもそも、ロバスト最適化とはロバスト性を最大にすることです。
では、ロバスト性を最大にするにはどうすればいいでしょうか？
それは、感度を最小化することです。
学生
感度の最小化？
もう、ついていけません….
大丈夫です。わかりやすく説明します。
では、質量meが誤差因子だった場合の「誤差因子に対する伝達関数G_jhの感度」はどのような式になるか考えてみてください。
これは、構造最適化で扱った感度と同じ形式になります。なので、以下になります。
$$
\frac{\partial G_{jh}}{\partial me}
$$
構造最適化の記事で下記のように説明しました。
感度が高ければ、質量mを付加することで伝達関数G_jhを大きく変化させることができます。一方、感度が低ければ、質量mを付加しても伝達関数G_jhは変化しません。
つまり、誤差因子meに対する伝達関数G_jhの感度が高いと、誤差因子meのばらつくことで、伝達関数G_jhが大きく変化してしまいます。
従って、ロバスト性の最大化問題は、感度の最小化問題となります。
（ここからが本題です。）
では、感度を最小化するにはどうすればいいでしょうか？
答えは、「感度の感度」を求めて、感度が最小になるように設計変数mを決定することです。文章だけでの説明ではしんどいので、例題として、下記例題を扱いましょう。
スポンサーリンク
○例題
図3のような4自由度のバネマス系でm2に力Fが作用します。製造ばらつきによってm2が±30%もばらついてしまいます。そこでロバスト性を向上させるために、質量mをm1～m4のどこの位置に付加してもいいとします。このときの質量付加前後の伝達関数G_4,2を求めなさい。ここで、質量はm1=m2=m3=m4=1 [kg] 、バネ定数はk1=k2=k3=k4=10^3、mは最大で3 [kg]とする。
図3 例題
上記の例題では”m1～m4の質量が変化しても良い”ということになるので、m1～m4が設計変数（制御因子）となります。また、m2が誤差因子です。
ちなみに、ロバスト最適化問題では感度を”1次感度”と呼び、感度の感度を”2次感度”と呼びます。
繰り返しになりますが、誤差因子m2に対する伝達関数G_jhの1次感度は以下です。
$$
1次感度： \frac{\partial G_{jh}}{\partial me}
$$
次に、2次感度について考えます。
2次感度は、”誤差因子m2に対する伝達関数G_jhの1次感度”を設計変数（ここではm3とします）で偏微分すればよいことになります。従って、下式になります。
$$
2次感度： \frac{\partial^2 G_{jh}}{\partial me \partial m_3}
$$
ただし、下の伝達関数合成法の式を用いると、mの1変数なので、1次感度を求めるときにm=m2(誤差因子)を代入して偏微分すると、1次感度はm2の変数となり、m2以外の変数で偏微分すると2次感度は0になってしまいます。
$
newG_{jh}=G_{jh}-\frac{G_{ji}G_{ih}}{G_{ii}-\frac{1}{ω^2 m}}
$$
なので、この問題では誤差は伝達関数Gに含まれており、付加質量mで2次感度を下げる問題として考えます。（初心者向けに問題を簡単にしたと思ってください。）
例題を用いたロバスト最適化の説明
では、下記の例題を解きながら詳しく説明したいと思います。
○例題
図3のような4自由度のバネマス系でm2に力Fが作用します。製造ばらつきによってm3が±30%もばらついてしまいます。そこでロバスト性を向上させるために、質量mをm1～m4のどこの位置に付加してもいいとします。このときの質量付加前後の伝達関数G_4,2を求めなさい。ここで、質量はm1=m2=m3=m4=1 [kg] 、バネ定数はk1=k2=k3=k4=10^3、mは最大で3 [kg]とする。なお、ロバスト性を向上させたい周波数は1.8Hzである。
1次感度は下記になります。
$$
\frac{\partial G_{4,2}}{\partial m}=ω^2 G_{4,i} G_{i,2}
$$
2次感度を導出します。
$$
\frac{\partial^2 newG_{jh}}{\partial m^2} = \frac{-2 G_{ji} G_{ii} G_{ih}}{ω^2(m G_{ii}-\frac{1}{ω^2} )^3}
$$
m→0に近いときなので、
$$
\frac{\partial^2 G_{jh}}{\partial m^2} = 2 ω^4 G_{ji} G_{ii} G_{ih}
$$
つまり、この例題では以下の2次感度を用いることになります。
$$
\frac{\partial^2 G_{4,2}}{\partial m^2} = 2 ω^4 G_{4,i} G_{i,i} G_{i,2}
$$
m3が±30%変動するとどうなるのか
そもそも、m3が±30%変動するとどうなるのでしょうか？MATLABのプログラムを組んで確認してみましょう。結果は図4のようになります。図4の黒線が理想状態（ばらつきのない状態）です。赤線がばらついたときの変動を表しています。青点線の位置が1.8Hzです。
図4 質量m3の変動と伝達関数の関係
m3がばらつくことで、伝達関数がばらつくことがわかります。
2次感度から1次感度を予測する
前回の復習になりますが、微小質量dmを付加したときの伝達関数newGは下式で予測できましたね。
$$
newG_{jh} =G_{jh} + \frac{\partial G_{jh}}{\partial m} dm
$$
この考え方と同様に2次感度を用いて1次感度を予測することができます。1次感度の予測式は以下になります。
$$
\frac{\partial newG_{jh}}{\partial m} =\frac{\partial G_{jh}}{\partial m} + \frac{\partial^2 G_{jh}}{\partial m^2} dm
$$
付加質量後の1次感度を最小化することがロバスト最適化です。
では1次感度を最小化したときのロバスト最適化結果を見てみましょう。
図5 ロバスト最適化結果
1.8Hzの1次感度を最小化する
↓
1次共振周波数を低下させる
↓
ロバスト最大化
という結果になりました。
「これでロバスト最適化と言われてもな～…..」という感じですよね。
私もそう思います。
今後、中～上級編のロバスト最適化について説明したいと思います。
もしリクエスト頂いて、その内容で記事が書けそうなら書いてみたい気もしています。
今日はこのへんでGood luck!
（ちなみに、この言葉はYoutuberのバフェット太郎さんのパクリです。）
スポンサーリンク
MATLABプログラム
コード
clear all
close all
% % % % 初期設定
k_vec=[1 1 1 1]*10^3;
[K]=eval_Kmatrix(k_vec);
freq=0.5:0.1:20;
f=[0;1;0;0];
C=zeros(length(k_vec));
% % % % m3がばらつく場合
dm3=(rand(1,100)-0.5)/0.5*0.3; %m3が30%ばらつくので、その影響を確認する
for ii1=1:length(dm3)
m_vec=[1 1 1+dm3(ii1) 1];
[M]=eval_Mmatrix(m_vec);
before_x=eval_direct_x(M,K,C,f,freq);
% % % % % %
figure(1)
loglog(freq,abs(before_x(4,:)),'r') %力fの成分が1なので、xを伝達関数Gは同じになるのでxを表示
hold on
end
% % % % 質量付加前の質量行列
m_vec=[1 1 1 1];
[M]=eval_Mmatrix(m_vec);
before_x=eval_direct_x(M,K,C,f,freq);
% % % % % %
figure(1)
loglog(freq,abs(before_x(4,:)),'k','linewidth',3) %力fの成分が1なので、xを伝達関数Gは同じになるのでxを表示
hold on
xlim([freq(1) freq(end)])
loglog([1.8 1.8],[10^-8 10^6],'b--')
xlabel('周波数 Hz')
ylabel('伝達関数 m/N')
ylim([10^-8 10^6])
% % % % 伝達関数の計算
G_2i_1p8Hz=eval_direct_x(M,K,C,f,1.8);
G_i4_1p8Hz=eval_direct_x(M,K,C,[0;0;0;1],1.8); %Gih 相反定理
% % % 自己応答Giiを求めるために計算
G_i3_1p8Hz=eval_direct_x(M,K,C,[0;0;1;0],1.8); %Gih 相反定理
G_i1_1p8Hz=eval_direct_x(M,K,C,[1;0;0;0],1.8); %Gih 相反定理
G_ii=[
G_i1_1p8Hz(1)
G_2i_1p8Hz(2)
G_i3_1p8Hz(3)
G_i4_1p8Hz(4)
];
% % % % 1次感度
dGdm=(2*pi*5)^2*G_2i_1p8Hz.*G_i4_1p8Hz;
% % % % 2次感度
d2Gdm2=2*(2*pi*1.8)^4*G_2i_1p8Hz.*G_i4_1p8Hz.*G_ii;
% % % % 2次感度から1次感度を予測する
% new_dGdm=dGdm + d2Gdm2*dm;
%
%
% % % % % % % 感度の算出
% G_2i_5Hz=eval_direct_x(M,K,C,f,5); %ややこしいので再計算
% G_i4_5Hz=eval_direct_x(M,K,C,[0;0;0;1],5); %Gih 相反定理
%
% dGdm=(2*pi*5)^2*G_2i_5Hz.*G_i4_5Hz;
%
% % % % % % % 最小化問題
% % m=[-0.9 2]
%
% % % % % % 質量付加後の質量行列
vec=[1 1 1 1];
dm=0.001; %微小な質量
target_f=1.8;
while sum(vec)<7
[tempdGdm,id]=sort(abs(d2Gdm2),'descend');
if id(1)==1
m_vec1=vec+[dm 0 0 0];
% m_vec2=vec+[-dm 0 0 0];
elseif id(1)==2
m_vec1=vec+[0 dm 0 0];
% m_vec2=vec+[0 -dm 0 0];
elseif id(1)==3
m_vec1=vec+[0 0 dm 0];
% m_vec2=vec+[0 0 -dm 0];
elseif id(1)==4
m_vec1=vec+[0 0 0 dm];
% m_vec2=vec+[0 0 0 -dm];
else
error('なにかおかしい')
end
[M1]=eval_Mmatrix(m_vec1);
% [M2]=eval_Mmatrix(m_vec2);
after_x1=eval_direct_x(M1,K,C,f,target_f);
% after_x2=eval_direct_x(M2,K,C,f,target_f);
if sum(m_vec1) >=7
break
end
vec=m_vec1;
M=M1;
% % % % % % 感度の再計算
G_2i_1p8Hz=eval_direct_x(M,K,C,f,1.8);
G_i3_1p8Hz=eval_direct_x(M,K,C,[0;0;0;1],1.8); %Gih 相反定理
% % % 自己応答Giiを求めるために計算
G_i3_1p8Hz=eval_direct_x(M,K,C,[0;0;1;0],1.8); %Gih 相反定理
G_i1_1p8Hz=eval_direct_x(M,K,C,[1;0;0;0],1.8); %Gih 相反定理
G_ii=[
G_i1_1p8Hz(1)
G_2i_1p8Hz(2)
G_i3_1p8Hz(3)
G_i4_1p8Hz(4)
];
% % % % 1次感度
dGdm=(2*pi*5)^2*G_2i_1p8Hz.*G_i4_1p8Hz;
% % % % 2次感度
d2Gdm2=2*(2*pi*1.8)^4*G_2i_1p8Hz.*G_i4_1p8Hz.*G_ii;
end
% % % % 最適化後の結果
% % % % % % % m3がばらつく場合
dm3=(rand(1,100)-0.5)/0.5*0.3; %m3が30%ばらつくので、その影響を確認する
for ii1=1:length(dm3)
m_vec= vec + [0 0 0+dm3(ii1) 0];
[M]=eval_Mmatrix(m_vec);
after_x=eval_direct_x(M,K,C,f,freq);
% % % % % %
figure(2)
loglog(freq,abs(after_x(4,:)),'r') %力fの成分が1なので、xを伝達関数Gは同じになるのでxを表示
hold on
end
% % % % 最適化後の結果
[M]=eval_Mmatrix(vec);
after_x=eval_direct_x(M,K,C,f,freq);
% % % % % %
figure(2)
loglog(freq,abs(after_x(4,:)),'k','linewidth',3) %力fの成分が1なので、xを伝達関数Gは同じになるのでxを表示
hold on
loglog([1.8 1.8],[10^-8 10^6],'b--')
ylim([10^-8 10^6])
xlim([freq(1) freq(end)])
xlabel('周波数 Hz')
ylabel('伝達関数 m/N')
○functionファイル
function [M]=eval_Mmatrix(m_vec) M=diag(m_vec); end
function [K]=eval_Kmatrix(k_vec)
K=zeros(length(k_vec));
for ii1=1:length(k_vec)
if ii1==1 K(ii1,ii1)=k_vec(ii1);
else
K(ii1-1:ii1,ii1-1:ii1)=K(ii1-1:ii1,ii1-1:ii1)+[1 -1;-1 1]*k_vec(ii1);
end
end
end
function x=eval_direct_x(M,K,C,f,freq)
x=zeros(length(f),length(freq));
j=sqrt(-1);
for ii1=1:1:length(freq)
w=2*pi*freq(ii1);
x(:,ii1)=inv(K+j*w*C-w^2*M)*f;
end
end
コメント