このページは機械翻訳を使用して翻訳されました。元の英語を参照するには、ここをクリックします。
カスタム並列シミュレーションのためのベクトル化サロゲート最適化
この例では、surrogateopt
UseVectorized
オプションを使用してカスタム並列最適化を実行する方法を示します。UseParallel
オプションを正常に使用できない場合に、この手法を使用できます。たとえば、UseParallel
オプションは、並列評価に parsim
を必要とする Simulink® シミュレーションには適用されない可能性があります。ベクトル化された並列シミュレーションを最適化するとかなりのオーバーヘッドが発生するため、この手法は時間のかかるシミュレーションに最も役立ちます。
この例の並列戦略は、最適化をサイズ N
のチャンクに分割することです。ここで、N
は並列ワーカーの数です。この例では、Simulink.SimulationInput
ベクトルに N
セットのパラメータを準備し、そのベクトルに対して parsim
を呼び出します。すべての N
シミュレーションが完了すると、 surrogateopt
はサロゲートを更新し、別の N
セットのパラメータを評価します。
モデルシステム
この例では、ローレンツ力学システムを短い時間間隔での等速円運動に適合させようとします。ローレンツ系とその均一円近似については、例 常微分方程式 (ODE) の適合 で説明します。
Lorenz_system.slx
Simulink モデルは Lorenz ODE システムを実装します。このモデルは、ライブ スクリプトを使用してこの例を実行するときに含まれます。
この例の最後にある fitlorenzfn
ヘルパー関数は、等速円運動からポイントを計算します。例 常微分方程式 (ODE) の適合 から、ローレンツ ダイナミクスにかなりよく一致する円運動パラメータを設定します。
x = zeros(8,1); x(1) = 1.2814; x(2) = -1.4930; x(3) = 24.9763; x(4) = 14.1870; x(5) = 0.0545; x(6:8) = [13.8061;1.5475;25.3616];
このシステムはシミュレーションにそれほど時間がかからないため、最適化の並列時間は、逐次最適化の時間よりも短くなりません。この例の目的は、並列でより適切に実行される特定の例を提供することではなく、ベクトル化された並列シミュレーションを作成する方法を示すことです。
目的関数
目的関数は、0 から 1/10 までの一連の時間にわたって、ローレンツ システムと等速円運動の差の二乗和を最小化することです。xdata
倍の場合、目的関数は
objective = sum((fitlorenzfn(x,xdata) - lorenz(xdata)).^2) - (F(1) + F(end))/2
ここで、lorenz(xdata)
は、時刻 xdata
におけるローレンツ システムの 3 次元進化を表し、F
は、円システムとローレンツ システムの対応する点間の距離の二乗のベクトルを表します。目的は、エンドポイントの値の半分を減算して、積分を最もよく近似することです。
一致する曲線として等速円運動を考慮し、シミュレーション内のローレンツパラメータを変更して目的関数を最小化します。
特定のパラメータのローレンツシステムを計算する
ローレンツの元のパラメータのローレンツ システムを計算してプロットします。
model = 'Lorenz_system'; open_system(model); in = Simulink.SimulationInput(model); % params [X0,Y0,Z0,Sigma,Beta,Rho] params = [10,20,10,10, 8/3, 28]; % The original parameters Sigma, Beta, Rho in = setparams(in,model,params); out = sim(in); yout = out.yout; h = figure; plot3(yout{1}.Values.Data,yout{2}.Values.Data,yout{3}.Values.Data,'bx'); view([-30 -70])
等速円運動を計算する
ローレンツ計算で時間間隔にわたって以前に与えられた x
パラメータの等速円運動を計算し、その結果をローレンツ プロットとともにプロットします。
tlist = yout{1}.Values.Time; M = fitlorenzfn(x,tlist); hold on plot3(M(:,1),M(:,2),M(:,3),'kx') hold off
この例の最後にある objfun
ヘルパー関数は、ローレンツ システムと等速円運動の差の二乗和を計算します。目的はこの二乗和を最小化することです。
ssq = objfun(in,params,M,model)
ssq = 26.9975
ローレンツシステムを並列に当てはめる
適合を最適化するには、surrogateopt
を使用して Simulink モデルのパラメータを変更します。この例の最後にある parobjfun
ヘルパー関数は、パラメーターのマトリックスを受け入れます。マトリックスの各行は、1 つのパラメーター セットを表します。この関数は、この例の最後にある setparams
ヘルパー関数を呼び出して、Simulink.SimulationInput
ベクトルのパラメータを設定します。次に、parobjfun
関数は parsim
を呼び出して、パラメータに基づいてモデルを並列に評価します。
並列プールを開き、プール内のワーカーの数として N
を指定します。
pool = gcp('nocreate'); % Check whether a parallel pool exists if isempty(pool) % If not, create one pool = parpool; end N = pool.NumWorkers
N = 6
BatchUpdateInterval
オプションを N
に設定し、 UseVectorized
オプションを true
に設定します。これらの設定により、surrogateopt
は一度に N
ポイントを目的関数に渡します。初期点を、以前に使用したパラメータに設定します。これは、均一な円運動にかなりよく適合するためです。MaxFunctionEvaluations
オプションを 600 に設定します。これは、この例で使用されているコンピューター上の 6 つのワーカーの整数倍です。
options = optimoptions('surrogateopt','BatchUpdateInterval',N,... 'UseVectorized',true,'MaxFunctionEvaluations',600,... 'InitialPoints',params);
現在のパラメータの上下 20% の範囲を設定します。
lb = 0.8*params; ub = 1.2*params;
速度を上げるには、シミュレーションを高速再起動を使用するように設定します。
set_param(model,'FastRestart','on');
目的関数の N
シミュレーション入力のベクトルを作成します。
simIn(1:N) = Simulink.SimulationInput(model);
再現性のためにランダムストリームを設定します。
rng(1)
parobjfun
を呼び出して、ベクトル化された並列方式で目的関数を最適化します。
tic [fittedparams,fval] = surrogateopt(@(params)parobjfun(simIn,params,M,model),lb,ub,options)
surrogateopt stopped because it exceeded the function evaluation limit set by 'options.MaxFunctionEvaluations'.
fittedparams = 1×6
10.3789 19.9803 10.0073 9.4956 2.8269 27.9263
fval = 23.0730
paralleltime = toc
paralleltime = 397.1877
目的関数の値が向上(減少)します。元の値と改善された値を表示します。
disp([ssq,fval])
26.9975 23.0730
適合した点をプロットします。
figure(h) hold on in = setparams(in,model,fittedparams); out = sim(in); yout = out.yout; plot3(yout{1}.Values.Data,yout{2}.Values.Data,yout{3}.Values.Data,'rx'); legend('Unfitted Lorenz','Uniform Motion','Fitted Lorenz') hold off
モデルを閉じるには、まず高速再起動を無効にする必要があります。
set_param(model,'FastRestart','off'); close_system(model)
まとめ
UseParallel
オプションを正常に使用できない場合は、surrogateopt
UseVectorized
オプションを true
に設定し、BatchUpdateInterval
オプションを並列ワーカーの数の倍数に設定することで、シミュレーションを並列で最適化できます。このプロセスは並列最適化を高速化しますが、オーバーヘッドが発生するため、時間のかかるシミュレーションに最適です。
補助関数
次のコードは補助関数 fitlorenzfn
を作成します。
function f = fitlorenzfn(x,xdata) theta = x(1:2); R = x(3); V = x(4); t0 = x(5); delta = x(6:8); f = zeros(length(xdata),3); f(:,3) = R*sin(theta(1))*sin(V*(xdata - t0)) + delta(3); f(:,1) = R*cos(V*(xdata - t0))*cos(theta(2)) ... - R*sin(V*(xdata - t0))*cos(theta(1))*sin(theta(2)) + delta(1); f(:,2) = R*sin(V*(xdata - t0))*cos(theta(1))*cos(theta(2)) ... - R*cos(V*(xdata - t0))*sin(theta(2)) + delta(2); end
次のコードは補助関数 objfun
を作成します。
function f = objfun(in,params,M,model) in = setparams(in,model,params); out = sim(in); yout = out.yout; vals = [yout{1}.Values.Data,yout{2}.Values.Data,yout{3}.Values.Data]; f = sum((M - vals).^2,2); f = sum(f) - (f(1) + f(end))/2; end
次のコードは補助関数 parobjfun
を作成します。
function f = parobjfun(simIn,params,M,model) N = size(params,1); f = zeros(N,1); for i = 1:N simIn(i) = setparams(simIn(i),model,params(i,:)); end simOut = parsim(simIn,'ShowProgress','off'); % Suppress output for i = 1:N yout = simOut(i).yout; vals = [yout{1}.Values.Data,yout{2}.Values.Data,yout{3}.Values.Data]; g = sum((M - vals).^2,2); f(i) = sum(g) - (g(1) + g(end))/2; end end
次のコードは補助関数 setparams
を作成します。
function pp = setparams(in,model,params) % parameters [X0,Y0,Z0,Sigma,Beta,Rho] pp = in.setVariable('X0',params(1),'Workspace',model); pp = pp.setVariable('Y0',params(2),'Workspace',model); pp = pp.setVariable('Z0',params(3),'Workspace',model); pp = pp.setVariable('Sigma',params(4),'Workspace',model); pp = pp.setVariable('Beta',params(5),'Workspace',model); pp = pp.setVariable('Rho',params(6),'Workspace',model); end
参考
surrogateopt
| parsim
(Simulink)