Optimization of resonant system based on Gaussian fit

Driver:

ActiveLearning

Download script: harmonic_oscillator_fit.m

The target of the study is to tune a resonant system. For simplicity, the resonant system is a harmonic oscillator with eigenfrequency \(\omega_0\) and damping \(\gamma\) driven with force \(F\) at frequency \(\omega\). The amplitude of the oscillator is

\[a_{\rm harm}(\omega) = \frac{F}{\sqrt{(2\omega\omega_0\gamma)^2 + (\omega_0^2 - \omega)^2}} .\]

We assume the system as a black box with unknown resonant behavior. To capture the resonant feature, it is fitted against a Gaussian plus a linear function

\[a_{\rm fit}(\omega) = A \exp\left(-\frac{1}{2}\frac{(\omega-\tau)^2}{\sigma^2}\right) + B\omega + C .\]

The target is to tune \(F, \omega_0, \gamma\) such that the fitted amplitude is \(A = 10\), the resonance frequency is \(\tau = 2\) and the resonance width \(\sigma\) is as small as possible.

jcm_optimizer_path = '<JCM_OPTIMIZER_PATH>';
addpath(fullfile(jcm_optimizer_path, 'interface', 'matlab'));
 
server = jcmoptimizer.Server();
client = jcmoptimizer.Client(server.port); 

%Amplitude of a driven harmonic oscillator. 
function a = amplitude(omega, F, omega0, gamma)
    a = F/sqrt((2*omega*omega0*gamma)^2 + (omega0^2 - omega^2)^2);
end

% Definition of the search domain to tune the oscillator
design_space = { ...
    struct('name', 'F', 'type', 'continuous', 'domain', [0.1, 40.0]), ... 
    struct('name', 'omega0', 'type', 'continuous', 'domain', [1.0, 4.0]), ...
    struct('name', 'gamma', 'type', 'continuous', 'domain', [0.03, 0.3]) ...
};

%The amplitude is scaned for 10 different driving frequencies
omegas = linspace(1,3,10);

%Since the system depends on the driving frequency which is not a design parameter,
%it is defined as an environment parameter.
environment = { ...
    struct('name', 'omega', 'type', 'variable', 'domain', [omegas(1), omegas(end)]) ... 
};

 % Creation of the study object with study_id 'harmonic_oscillator_fit'
study = client.create_study( ...
    'design_space', design_space, ...
    'environment', environment, ...
    'driver','ActiveLearning',...
    'name','Optimization of resonant system based on Gaussian fit',...
    'study_id', 'harmonic_oscillator_fit');

%definition of variables

%The single-output Gaussian process is scanned over all omegas
scan = struct();
scan.type = "Scan";
scan.name = "amplitude_scan";
scan.output_dim = length(omegas);
scan.scan_parameters = {"omega"};
scan.scan_values = omegas.';
scan.input_surrogate = "amplitude";

%The result of the omega-scan is fitted to a Gaussian + linear expression with
%amplitude A, resonance frequency tau, linear gradient B, and constant offset C.
%The output are the fitted values and the mean-squared error (MSE).
fit = struct();
fit.type = "Fit";
fit.name = "fit";
fit.input = "amplitude_scan";
fit.expression = "A*exp(-0.5*(omega-tau)^2/sigma^2) + B*omega + C";
fit.output_names = {"A", "B", "C", "tau", "sigma", "MSE"};
fit.output_dim = 6;
fit.model_variables = {"omega"};
fit.variable_values = omegas.';
fit.initial_parameters = [1.0, 0.0, 0.0, 1.0, 0.5];

%The expression that defines the loss function:
expression = struct();
expression.type = "Expression";
expression.name = "loss";
expression.expression = "(A-20)^2 + (tau - 2)^2 + sigma^2 + 0.001*MSE";
 
%configuration of study
study.configure( ...
    'max_iter', 30, ...
    'surrogates', { ...
        ...A single-output Gaussian process learns the dependence of the amplitude on the
        ...design and environment parameters.
        struct('type', "GP", 'name', "amplitude", 'output_dim', 1, 'tolerance', 1e-12) ...
    }, ...
    'variables', {scan, fit, expression}, ...
    'objectives', { ...
        ...The only objective of the study is to minimize the loss.
        struct("type", "Minimizer", "variable", "loss"), ...
    }, ...
    'acquisition_optimizer', struct( ...
        ...Sample computation based on fits is more expensive. In this case
        ...advanced sample computation is usually too laborious.
        'compute_suggestion_in_advance',  false ...
    ) ...
);

% Evaluation of the black-box function for specified design parameters
function observation = evaluate(study, sample)

    pause(2); % make objective expensive
    observation = study.new_observation();
    %The harmonic-oscillator amplitude is evaluated for all omega values and the
    %observed values are used as training input for the Gaussian process "amplitude"
    for i = 1:length(omegas)
        observation.add(amplitude(omegas(i), sample.F, sample.omega0, sample.gamma), ...
                        'environment_value', {omegas(i)}, 'model_name', "amplitude");
    end
    
end  

% Run the minimization
study.set_evaluator(@evaluate);
study.run(); 


% Print result
best = study.get_state("driver.best_sample");
fprintf("Best design parameters: F=%f, omega0=%f, gamma=%f\n", ...
        best.F, best.omega0, best.gamma);
fprintf(['Objective: (A-20)^2 + (tau - 2)^2 + sigma^2 + 0.001*MSE = ', ...
         '%.3f\n'], study.get_state('driver.min_objective'));