simulationサンプルコード¶
例1：初期値応答¶
- 解析内容
- 解析対象モデル：IEEE68busモデル
- シミュレーション時間： 0~20秒
- シミュレーション条件
母線2に接続された機器の3番目の状態が平衡点から0.1ズレた状態を初期値とした場合。IEEE68busモデルの場合母線3に接続されているのは
generator_1axisであり、このモデルの3番目の状態は周波数偏差に対応しているため、この条件は母線２に接続された同期発電機の周波数が何らかの影響でズレた場合の応答を調べるということになります。
% ネットワークの定義
net = network_IEEE68bus();
% シミュレーションのためのオプションを定義・決定
option = struct();
option.x0_sys = net.x_equilibrium;
option.x0_sys(16)= option.x0_sys(16) + 0.01;
% シミュレーションの実行
out = net.simulate([0 20], option);
net = network_IEEE68bus();
x0 = net.x_equilibrium;
x0(16)= x0(16) + 0.01;
out = net.simulate([0 20],'x0_sys',x0);
ちなみに、母線３に接続された機器の２番目の状態が、上のコードのように系統全体の状態の16番目に対応しているのは、系統全体の状態を母線番号順に機器の状態を繋げたベクトルとしておりIEEE68busモデルの場合、母線1と2に状態を7つ持つ機器が接続されているためです。しかし、系統内の何番目に状態を何個持つ機器が接続されているか把握するのは難しいでしょう。そういった場合は以下の様に各componentクラスのget_nxという関数を上手く使うと、より一般的な実装とすることが可能です。
bus_idx = 3; %3番目の母線
state_idx = 2; %2番目の状態
%母線1から母線(bus_idx-1)番目までの状態の個数を足し合わせる
idx = 0;
for i = 1:bus_idx-1
idx = idx+net.a_bus{i}.component.get_nx;
end
%状態の2番目のインデックスを知りたいため
idx = idx+state_idx;
%IEEE68busモデルの場合は、idx = 7+7+2= 16となります。
option.x0_sys(idx)= option.x0sys(idx) + 0.1;
：
例2：外乱応答¶
- 解析内容
- 解析対象モデル：IEEE68busモデル
- シミュレーション時間： 0~30秒
- シミュレーション条件 母線1に0〜0.07秒の間と15〜15.05秒の間に、母線2と母線5〜7に10~10.01秒の間、地絡が起きた場合
% ネットワークの定義
net = network_IEEE68bus();
% シミュレーションのためのオプションを定義・決定
option = struct();
option.fault = {{[0 0.07], 1}, {[15 15.05], 1}, {[10,10.01],[2,5:7]}};
% シミュレーションの実行
out = net.simulate([0 30], option);
net = network_IEEE68bus();
out = net.simulate([0 30],...
'fault',{{[0 0.07], 1}, {[15 15.05], 1}, {[10,10.01],[2,5:7]}});
例3：入力応答¶
(一連の解析実行例の解説ページと同一の解析です。)
- 解析内容
- 解析対象モデル：一連の解析実行例の解説ページで作成した3busモデル
- シミュレーション時間： 0~30秒
- シミュレーション条件：母線３に接続された負荷のインピーダンスが大きくなった場合。
% ネットワークの定義
net = network_sample3bus();
%条件設定
time = [0,10,20,60];
u_idx = 3;
u = [0, 0.05, 0.1, 0.1;...
0, 0, 0, 0];
%解析実行
out = net.simulate(time,u, u_idx);
例4：初期値・外乱・入出力応答（組み合わせ）¶
ここまでの設定を全て組み合わせて以下のような場合のシュミレーションを実行させます。「バス1の初期状態のみを平衡点から少しずらし，00.07秒にバス1，1010.05秒にバス10の地絡を外乱として印加する状況で，バス1,10にランダムの入力を加えた時の応答」
% ネットワークの定義
net = network_IEEE68bus();
% シミュレーションのためのオプションを定義・決定
option = struct();
option.x0_sys = net.x_equilibrium;
option.x0_sys(16)= option.x0_sys(16) + 0.01;
option.fault = {{[0 0.07], 1}, {[15 15.05], 1}, {[10,10.01],[2,5:7]}};
%条件設定
time = [0,10,20,60];
u_idx = 20;
u = [0, 0.05, 0.1, 0.1;...
0, 0, 0, 0];
% シミュレーションの実行
out = net.simulate(time, u, u_idx, option);
net = network_IEEE68bus();
x0 = net.x_equilibrium;
x0(16)= x0(16) + 0.01;
out = net.simulate([0 20],u,u_idx,...
'x0_sys',x0,...
'fault',{{[0 0.07], 1}, {[15 15.05], 1}, {[10,10.01],[2,5:7]}});
シミュレーション結果を見る¶
今回解析対象とした、IEEE68busモデルの発電機母線は1〜16番目となっています。また、これらの発電機母線の発電機モデルは全てgenerator_1axisクラスのインスタンスとなっています。そのため例えば各発電機院周波数偏差が見たいとなった場合は、各機器の2番目の状態を見れば良いということになります。このような場合は簡単で、以下の様にプロットすれば良いです。
figure;
hold on
arrayfun(@(idx) plot(out.t,out.X{idx}(:,2)),1:16);
xlabel('時刻(s)','FontSize',15);
ylabel('周波数偏差','FontSize',15);
legend()
title('各同期発電機の周波数偏差','FontSize',20)
hold off
simulationResultクラスを用いたプロットでは内部のmethodで各機器の状態変数名を用いて自動で分類を行いプロットします。
Tips¶
前ページでも取り扱いましたがGUILDAにはシミュレーション結果を可視化するためのクラスsimulationResult classが実装されています。上のsimulate実行の際の引数に'tools',trueを付け足すことで、出力結果をこのクラスのインスタンスへと変更されます。
- 例４の実行例でsiimulationResultクラスの出力にする場合
net = network_IEEE68bus();
x0 = net.x_equilibrium;
x0(16)= x0(16) + 0.1;
out = net.simulate([0 20],u,u_idx,...
'x0_sys',x0,...
'fault',{{[0 0.07], 1}, {[15 15.05], 1}, {[10,10.01],[2,5:7]}},...
'tools',true);
プロットを行う際は、以下の様に実行します。
out.UIplot