Optimal control of a system in a changing environment

Driver:

ActiveLearning

Download script: changing_environment.m

The target of the study is show how to control a system that depends on environment parameters such as temperature or humidity. While the environment parameters can be measured, their influence on the system’s performance is often unknown.

As an example objective the 2d Rastrigin function

\[r(x_1, x_2) = 10\cdot 2 +x_1^2 + x_2^2 - 10\cos(2\pi x_1) - 10\cos(2\pi x_1)\]

is considered. The environment parameter \(\phi\) acts as an additional phase offset to the first cosine function in the objective function

\[r(x_1, x_2, \phi) = 10\cdot 2 +x_1^2 + x_2^2 - 10\cos(2\pi x_1 + \phi) - 10\cos(2\pi x_1)\]

The phase shall slowly vary over time as

\[\phi(t) = 2\pi\sin\left(\frac{t}{3{\rm min}}\right).\]

Please note, that this specific time dependent behaviour is not exploited and is assumed to be unknown.

Before being able to control the system in an optimal way depending on the environment, one has to learn for many environment values, where the global minimum is located. To this end, a standard Bayesian optimization is performed for 500 iterations that explores the parameter space. In a second phase, the target is to evaluate the system in an optimal way, i.e. an exploration of the parameter space is not desired. This behaviour is mainly enforced by choosing a small scaling value.

The control phase could have an arbitrary number of iterations and it would be problematic to add all new observations to the study. On the one hand, this slows down the computation time of a suggestion. Since the environment value changes during the computation, this can lead to less optimal evaluation points. On the other had, adding more and more data points close to each other leads to an ill conditioned Gaussian process surrogate. To avoid these drawbacks, data points are not added in the control phase if the study predicts a value with very small uncertainty, which means that the observation would not add significant information.

client = jcmoptimizer.Client(); 

%Rastrigin-like function depending on additional phase offset phi 
function val=rast(x1, x2, phi)
    val = 10*2 + x1^2 + x2^2 - 10*cos(2*pi*x1 + phi) - 10*cos(2*pi*x2);
end
 
%time-dependent slowly varying phi
function phi=current_phi()
    t = now()*1e5;
    phi = 2*pi*sin(t/180);
end
 
% Definition of the search domain
design_space = { ...
    struct('name', 'x1', 'type', 'continuous', 'domain', [-1.5,1.5]), ...
    struct('name', 'x2', 'type', 'continuous', 'domain', [-1.5,1.5]) ...
};

% Definition of the environment variable "phi"
environment = { ...
    struct('name', 'phi', 'type', 'variable', 'domain', [-2*pi, 2*pi]) ...
};

 % Creation of the study object with study_id 'changing_environment'
study = client.create_study( ...
    'design_space', design_space, ...
    'environment', environment, ...
    'driver','ActiveLearning',...
    'study_name','Optimal control of a system in a changing environment',...
    'study_id', 'changing_environment');

% In the initial training phase, the target is to explore the
% parameter space to find the global minimim.
study.configure( ...
    ... train with 500 data points
    'max_iter', 500, ...    
    ... Advanced sample computation is switched off since the environment
    ... parameter phi can change significantly between computation
    ... of the suggestion and evaluation of the objective function
    'acquisition_optimizer', struct('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();
    % get current phi
    phi = current_phi();
    observation.add(rast(sample.x1, sample.x2, phi), 'environment_value', {phi});
    
end  

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


% The target in the control phase is to evaluate the offet Rastrigin function only
% at well performing (x1,x2)-point depending on the current value of the environment.
MAX_ITER = 500; %evaluate for 500 additional iterations
study.configure( ...
    'max_iter', 500 + MAX_ITER, ...
    ... The scaling is reduced to penalize parameters with large uncertainty    
    'scaling', 0.01, ...
    ... The lower-confidence bound (LCB) strategy is chosen instead of the
    ... default expected improvement (EI). LCB is easier to maximize at the
    ... risk of less exploration of the parameter space, which is anyhow not
    ... desired in the control phase.
    'objectives', {...
        struct('type', 'Minimizer', 'name', 'objective', 'strategy', 'LCB') ...
    }, ...
    'acquisition_optimizer', struct('compute_suggestion_in_advance', false) ...
);

% keep track of suggested design points and phis at request time and evaluation time
design_points = {};
phis_at_request = {};
phis_at_eval = {};
            
iter = 0;
while not(study.is_done())
    iter = iter + 1;
    if iter > MAX_ITER
        break
    end
        
    phi = current_phi();
    suggestion = study.get_suggestion({phi});
    phis_at_request{end + 1} = phi;
    sample = suggestion.sample;
    design_points{end + 1} = [sample.x1, sample.x2]; 
    try
        obs = evaluate(study, sample);
        %update phi from observation
        phi = obs.data.None{1}.env{1};
        phis_at_eval{end + 1} = phi;
        
        predictions = study.driver.predict({[sample.x1, sample.x2, phi]});
        std = sqrt(predictions.variance(1, 1));
            
        fprintf("Uncertainty of prediction %f\n", std);    
        % add data only if prediction has significant uncertainty
        if std > 0.01
            study.add_observation(obs, suggestion.id);
        else
            study.clear_suggestion( ...
                suggestion.id, sprintf("Ignoring observation with uncertainty %f", std));
        end
            
    catch err
        study.clear_suggestion( ...
            suggestion.id, ...
            sprintf("Evaluator function failed with error: %s", err.message) ...
        )
        rethrow(err);
    end    

end
            

fig = figure('Position', [0, 0, 1000, 500]);
 
% all observed training samples
observed = study.driver.get_observed_values();
subplot(1, 2, 1)
plot(observed.means,".")
hold on;
xline(500, '--');
xlabel("training+control iteration")
ylabel("observed value of Rastrigin function")

%observed values during control phase
observed_vals = [];

%values that would have been observed at request time,
%i.e. if there would be no time delay between request and 
%evaluation of suggestion
observed_vals_at_request = [];

for i = 1:length(design_points)
    observed_vals(i) = rast( ...
        design_points{i}(1), design_points{i}(2), phis_at_eval{i});
    observed_vals_at_request(i) = rast( ...
        design_points{i}(1), design_points{i}(2), phis_at_request{i});
end

%best value of x1-parameter depending on environment
function val = best_x1(phi)
    val = -phi/(2*pi);
    if abs(phi) > pi
        val = val + sign(phi);
    end
end

% best possible values 
best_vals = [];
for i = 1:length(phis_at_eval)
    phi = phis_at_eval{i};
    best_vals(i) = rast(best_x1(phi), 0, phi);
end

subplot(1, 2, 2)
semilogy(observed_vals,".", 'DisplayName', 'observed values')
hold on;
semilogy(observed_vals_at_request,".", 'DisplayName', 'observed values if no time delay')
semilogy(best_vals, 'DisplayName', 'best possible values')
ylim([1e-4, 1e1])
xlabel("control iteration")
legend()
saveas(fig, "training_and_control.png")

client.shutdown_server();

Left: During the initial training phase in the first 500 iterations, the parameter space is explored leading to small and large objective values. In the control phase, only small objective values are observed. Right: The observed values (blue dots) agree well with the lowest achievable values (green line). Most of the deviations are due to the time offset between the request of a new suggestion for a given environment value \(\phi\) and the actual evaluation of the Rastrigin function about a second later. To see this, the values that would have been observed at the time of request are shown as orange dots.