Модельная система

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

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

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

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 параметр to 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