Система модели

Этот пример пытается соответствовать Лоренцу динамическая система к универсальному круговому движению на кратковременном интервале. Система Лоренца и ее универсальное круговое приближение описаны в Подгонке в качестве примера Обыкновенное дифференциальное уравнение (ODE).

Lorenz_system.slx Модель Simulink реализует систему ОДУ Лоренца. Эта модель включена, когда вы запускаете этот пример с помощью live скрипта.

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) представляет 3-D эволюцию системы Лоренца во времена xdata, и 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 функция помощника в конце этого примера принимает матрицу параметров, где каждая строка матрицы представляет один набор параметров. Вызовы функции 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

Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 6).

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(100)

Оптимизируйте цель векторизованным параллельным способом путем вызова 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.5627   19.8962    9.8420    8.9616    2.5723   27.9687

fval = 23.6361

paralleltime = toc

paralleltime = 457.9271

Значение целевой функции улучшается (уменьшается). Отобразите исходные и улучшенные значения.

disp([ssq,fval])

   26.9975   23.6361

Постройте подходящие точки.

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