このページは機械翻訳を使用して翻訳されました。元の英語を参照するには、ここをクリックします。

カスタム並列シミュレーションのためのベクトル化サロゲート最適化

この例では、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)

関連するトピック

• 代理最適化