Моделируйте систему

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

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

The fitlorenzfn вспомогательная функция в конце этого примера вычисляет точки из равномерного кругового движения. Установите параметры кругового движения из примера Подгонка обыкновенного дифференциального уравнения (ОДУ), которые достаточно хорошо соответствуют динамике Лоренца.

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

The objfun вспомогательная функция в конце этого примера вычисляет сумму различий между системой Лоренца и равномерным круговым движением. Цель состоит в том, чтобы минимизировать эту сумму квадратов.

ssq = objfun(in,params,M,model)

ssq = 26.9975

Подгонка системы Лоренца параллельно

Чтобы оптимизировать подгонку, используйте surrogateopt для изменения параметров модели Simulink. The parobjfun вспомогательная функция в конце этого примера принимает матрицу параметров, где каждая строка матрицы представляет один набор параметров. Функция вызывает setparams вспомогательная функция в конце этого примера, чтобы задать параметры для Simulink.SimulationInput вектор. The 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