Contents

Neuromechanic - A neuromechanic model of button-pressing

Uses: bayesopt, Simulink, Simscape

Authors: Antti Oulasvirta (1), Byungjoo Lee (2), Sunjun Kim (1) (1) Aalto University; (2) KAIST

Version: September 19 2017

This is the version reported in the Proceedings of CHI 2018 (Oulasvirta, Kim, & Lee, 2018). Read full paper at http://doi.org/10.1145/3173574.3174082

This work is licensed under the Creative Commons Attribution 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.

warning('off','all');   % All warnings off
close all;              % Clears workspace
addpath(genpath(pwd));  % Adds working folder to path

Simulink and Simscape initialization

CFR_HomeDir = pwd;
CFL_libname = 'Libraries/Contact_Forces_Lib';
load_system(CFL_libname);
load_system('CMFB');

% "hws" means "human workspace", used for storing variables in the human model
hws = get_param('CMFB', 'modelworkspace');
hws.DataSource = 'MAT-File';
hws.FileName = 'params';
hws.assignin('peak_int', 0.5);  % Set peak of integrated p-center at 500 ms (this will be udpated by the BO with experience)
hws.saveToSource;
hws.reload;

Bayesian Optimizer setup

See Matlab documentation for reference

max_Time = 1200;                                          % Optimization budget (in seconds)
BOexplorationratio = 0.15;                              % Exploration rate (larger value -> more exploration)
traintrials = 8;                                        % Number of repeated attempts of a button press
testtrials = 100;                                       % Number of test trials (after model-learning)

Motor, finger and sensor feedback model parameters

See Table 1 in paper for explanation

Motornoise = readtable('i_motornoise.csv');             % Muscle noise parameters
Sensory = readtable('i_sensorynoise.csv');              % Sensory feedback noise parameters
Finger = readtable('i_finger.csv');                     % Physical properties of finger
Finger_Stiffness = readtable('i_finger_stiffness.csv'); % Finger pulp model (stiffness per displacement)
Finger.flength = 0.06;                                  % Length of finger bone (unit: meters)

% Hill-type finger muscle model parameters
Finger.Lmag = Finger.L0ag * 0.3;                        % Minimum length at which the muscle can produce force (agonist)
Finger.Lman =Finger.L0an * 0.3;                         % Minimum length at which the muscle can produce force (antagonist)
Finger.Ks = 0.8 * Finger.fcsa;                          % Stiffness coefficient of passive muscle element
Finger.Kd = 0.1 * Finger.Ks;                            % Damping coefficient of passive muscle element
Finger.vmag = Finger.L0ag * 8;                          % Maximum contraction velocity under no load (agonist)
Finger.vman = Finger.L0an * 8;                          % Maximum contraction velocity under no load (antagonist)

Button model parameters

Button = readtable('i_button.csv');             % Muscle noise parameters
Feedback=readtable('i_feedback.csv');                   % Physical properties - intensity, loudness - of of light and audio feedback from button-press, respectively
Button_Stiffness_temp = readtable('FDdatabase/Mechanical/Samsung SKM-1000UB (2015-05, Cherry MX Blue).csv'); % Tactile button type
%Button_Stiffness_temp = readtable('o_linear_fd.csv');  % Linear button
d = designfilt('lowpassfir', ...
    'PassbandFrequency',0.15,'StopbandFrequency',0.2, ...
    'PassbandRipple',1,'StopbandAttenuation',60, ...
    'DesignMethod','equiripple');
%Button_Stiffness_temp.Var2 = filtfilt(d,Button_Stiffness_temp.Var2); %Smoothing the empirical stiffness of buttons; Comment this out for Samsung tactile button
maximum_travel_index=find(Button_Stiffness_temp.Var1==max(Button_Stiffness_temp.Var1));
x_displacement=[Button_Stiffness_temp.Var1(1:maximum_travel_index)/1000000;(max(Button_Stiffness_temp.Var1)/1000000+0.0001)];
stiffness=[0;diff(Button_Stiffness_temp.Var2(1:maximum_travel_index)/100)./diff(Button_Stiffness_temp.Var1(1:maximum_travel_index)/1000000);10000];
activation_d=zeros(size(stiffness));
activation_d(1)=-0.002;                                 % Activation point of the physical button
activation_midair=zeros(size(stiffness));
activation_midair(1)=-100000;                           % Disable the activation of mid-air button
activation_p=zeros(size(stiffness));
activation_p(1)=100000;                                 % Disable the activation of touch button
damping=zeros(size(stiffness));
damping(1)=0.1;                                         % Damping of the button spring
Button_Stiffness=table(x_displacement,stiffness,activation_p,activation_d,activation_midair,damping);
Finger.fstart_height = Finger.fradius+Button.bdepth/2;  % Starting height of finger
Finger.fstart_height_std = 0.001;                       % Standard deviation (uncontrollable component) of finger starting height

OptimizableVariables for Bayesian Optimization

These variables define the decision task of the BO, ranges in square brackets Motor command: $\theta = \{ \mu_{A+}, \tau_{A+}, \sigma_{A+}\}$ % Expected p-center: $ pc_e $ % The current setup assumes a 3.0 second total time window for the simulation

onag = optimizableVariable('onag', [0.0, 1.5]);         % Onset of agonist muscle (s)
amag = optimizableVariable('amag', [0.5, 1.0]);         % Amplitude of agonist
pwag = optimizableVariable('pwag', [1, 5]);             % Pulse width of agonist
exp_ap = optimizableVariable('exp_ap', [0, (3-0.1)]);   % Expected p-center

Training phase

Obtains an optimized motor program over repeated button presses Uses Expected Improvement as the Acquisition function of BO See MATLAB's documentation for description of BO's parameters

disp('Training');
testing = false;
training = @(x)perceptualcontrol(x,traintrials,Motornoise,Finger,Feedback, Button,Sensory, testing);
results = bayesopt(training, [amag,pwag,onag,exp_ap],'MaxTime',max_Time, 'Verbose', 1, 'PlotFcn',[], ...
    'ExplorationRatio', BOexplorationratio, 'NumSeedPoints', 5,'IsObjectiveDeterministic',false, ....
    'MaxObjectiveEvaluations', inf, ...
    'AcquisitionFunctionName', 'expected-improvement');
Training

Testing phase

Uses the optimized motor program over repeated trials for estimates of performance Results are written in /results

writetable (results.XAtMinObjective,'results/trace');   % Stores winner motor program
disp('Testing')
testing = @(x,testtrials)perceptualcontrol(x,testtrials,Motornoise,Finger,Feedback,Button,Sensory, true); % Testing = true
objvalue = testing(results.XAtMinObjective,testtrials);
Testing

Perceptual control system linked to Simulink and Simscape

This system runs the given motor program and returns objective value to BO See CMFB.slx for the finger + button models

function [objective] = perceptualcontrol(x,n,Neural,Finger,Feedback,Button,Sensory,testing)
    % Load the hws Simulink model, access variables in modelworkspace
    hws = get_param('CMFB', 'modelworkspace');
    hws.DataSource = 'MAT-File';
    hws.FileName = 'params';

    % Passive record-keeping of p-centers for final reporting
    peak_ps = [];
    peak_ds = [];
    peak_ls = [];
    peak_as = [];

    % Passive record-keeping of activation and other data for final
    % reporting
    activation_points = [];                                 % Button activations
    activations = [];                                       % Vector of activations (button not activated)
    ovs = [];                                               % Vector of objective values obtained in repeated trials
    deltas = [];                                            % Vector of deltas between expected and obtained perceptions
    percepts = [];                                          % Vector of percepts (after button presssing)
    displacements = [];
    velocities = [];
    muscleforces = [];

    % n Repeated button-pressing trials, over which the objective score is
    % computed
    for i = 1:n

        % Muscle model: Set trial-specific input parameters

        % Noise model: additive Gaussian noise
        hws.assignin('amplitude_ag', abs(random('norm',x.amag,Neural.ag_am_noise)));    % Amplitude of agonist
        hws.assignin('pulsewidth_ag', abs(random('norm',x.pwag,Neural.ag_pw_noise)));   % Pulse width of agonist
        hws.assignin('amplitude_an', 0);
        hws.assignin('pulsewidth_an', 10);                  % In the CHI2018 version, antagonist is passive (resists but not activated)
        ag_pd = min(3,random('norm',x.onag+Neural.ag_motor_delay,Neural.ag_motor_delay_noise));
        hws.assignin('phasedelay_ag', ag_pd);               % Onset time is subjected to noise and delay caused by to Treisman gating
        finger_start_height = random('norm',Finger.fstart_height,Finger.fstart_height_std); % Varying start position of finger
        hws.assignin('finger_start_height', finger_start_height);
        hws.assignin('window', 3);                          % Sets the total time window of the simulation (3 seconds in paper)
        hws.saveToSource;
        hws.reload;                                         % Save and refresh

        % Simulate
        sim('CMFB');


        % The rest of this function is used for computing objective score
        % and reporting variables on physical events in the simulation

        % Muscle force
        muscleforce = abs(min(push_muscle_force.data)); % amount of force
        muscleforces = [muscleforces muscleforce];

        % Activation point
        max_ap = max(keyevent.data);
        peak_ap = find(keyevent.data == max_ap);  % the first keyevent
        peak_ap = keyevent.Time(peak_ap(1));

        size_d = size(raw_displacement.data);
        size_v = size(velocity.data);

        % Computing p-centers
        % Pressure sensation in finger pulp
        % Note that sensory noises are applied in CMBF.slx
        max_p = max(pulpcontact.data);              % amplitude of peak
        peak_p = find(pulpcontact.data == max_p);   % location of first peak
        peak_index=peak_p(1);                       % find first peak
        peak_p = pulpcontact.Time(peak_index);      % take the value at that peak
        peak_p_shifted = peak_p + Sensory.p_delay;  % shift the peal per sensory delay
        peak_p_asdn = peak_p_shifted;
        p_power=max_p;
        p_reliability=1/log2(p_power/Sensory.p_noise);
        peak_ps = [peak_ps peak_p_asdn];

        % Percept: Fingertip displacement
        min_d = min(displacement.data);
        peak_d = find(displacement.data == min_d);
        peak_index=peak_d(1);
        peak_d = displacement.Time(peak_index);
        peak_d_shifted = peak_d + Sensory.d_delay;
        peak_d_asdn = peak_d_shifted;
        peak_ds = [peak_ds peak_d_asdn];
        d_power=Finger.fstart_height-min_d;
        d_reliability=1/log2(d_power/Sensory.d_noise);

        % Percept: Light
        max_l = max(keyeventtime_l.data);
        peak_l = find(keyeventtime_l.data == max_l);  % the first light
        peak_l = keyeventtime_l.Time(peak_l(1));
        peak_l_shifted = peak_l +Sensory.l_delay;
        peak_l_asdn = peak_l_shifted;
        peak_ls = [peak_ls peak_l_asdn];
        l_power=Feedback.l_intensity*max_l;
        l_reliability=1/log2(l_power/Sensory.l_noise);

        % Percept: Auditory
        max_a = max(keyeventtime_a.data);
        peak_a = find(keyeventtime_a.data == max_a);  % the first beep
        peak_a = keyeventtime_a.Time(peak_a(1));
        peak_a_shifted = peak_a +Sensory.a_delay;
        peak_a_asdn = peak_a_shifted;
        peak_as = [peak_as peak_a_asdn];
        a_power=Feedback.a_loudness*max_a;
        a_reliability=1/log2(a_power/Sensory.a_noise);

        n_keyevents = sum(keyevent.data);

        % Objective function is set to 0 and increased based on success and
        % precision and effort
        objectivevalue = 0.0;

        raw_max_p = max(raw_pulpcontact.data);

        % Cue integration MLE
        % Because active sensory inputs are different per button type, we
        % separate some parts of computation
        if(n_keyevents==1&&raw_max_p~=0) % Physical and touch button single activation
            sum_var = 1/p_reliability + 1/d_reliability + 1/l_reliability + 1/a_reliability;
            peak_int = (1 / p_reliability / sum_var) * peak_p_asdn + (1 / d_reliability / sum_var) * peak_d_asdn + (1 / l_reliability / sum_var) * peak_l_asdn + ((1 / a_reliability) / sum_var) * peak_a_asdn;
            activations = [activations 1];
            activation_points = [activation_points peak_ap];
        elseif (n_keyevents==0&&raw_max_p~=0) % Physical button contacted but not activated
            sum_var = 1/p_reliability + 1/d_reliability;
            peak_int = ((1 / p_reliability) / sum_var) * peak_p_asdn + ((1 / (d_reliability)) / sum_var) * peak_d_asdn;
            objectivevalue = objectivevalue + 1; % + 0.5 / muscleforce;
            % HACK: encourages the use of more force to activate the button
            activations = [activations 0];
        elseif (n_keyevents==1&&raw_max_p==0) % Midair button not contacted but activated
            sum_var = 1/d_reliability + 1/l_reliability + 1/a_reliability;
            peak_int = (1 / d_reliability / sum_var) * peak_d_asdn + (1 / l_reliability / sum_var) * peak_l_asdn + ((1 / a_reliability) / sum_var) * peak_a_asdn;
            activations = [activations 1];
            activation_points = [activation_points peak_ap];
        else % If not contacted and not activated
            sum_var = 1/d_reliability;
            peak_int = ((1 / (d_reliability)) / sum_var) * peak_d_asdn;
            objectivevalue = objectivevalue + 2; % + 0.5 / muscleforce;
            activations = [activations 0];
        end

        percepts = [percepts peak_int];
        delta = abs(peak_int - x.exp_ap);
        deltas = [deltas delta];
        objectivevalue = objectivevalue + delta;
        ovs = [ovs objectivevalue];

    end

Print result plots to folder /results

    if testing == true % Reporting

        delta_avg = mean(deltas);
        delta_std = std(deltas);
        success_avg = mean(activations);
        activation_point_std = std(activation_points);
        activation_point_avg = mean(activation_points);
        ovs_avg = mean(ovs);
        muscleforce_avg = mean(muscleforces);
        muscleforce_std = std(muscleforces);
        results = table(delta_avg,delta_std,success_avg,ovs_avg,activation_point_avg,activation_point_std,muscleforce_avg,muscleforce_std);
        disp(results);
        writetable(results,'results/results.txt');

        plot(pulpcontact.time, pulpcontact.data,'Color','black');
        title('Pulp contraction');
        saveas(gcf,'results/pulpcontraction.fig')
        saveas(gcf,'results/pulpcontraction.png')

        plot(displacement.time, displacement.data,'Color','black');
        title('Finger tip position');
        saveas(gcf,'results/tipposition.fig');
        saveas(gcf,'results/tipposition.png');

        plot(velocity.time, velocity.data,'Color','black');
        title('Finger tip velocity');
        saveas(gcf,'results/tipvelocity.fig');
        saveas(gcf,'results/tipvelocity.png');

        plot(push_muscle_force.time, push_muscle_force.data,'Color','black');
        title('Muscle force');
        saveas(gcf,'results/muscleforce.fig');
        saveas(gcf,'results/muscleforce.png');

        plot(push_muscle_force.time, abs(push_muscle_force.data),'Color','black');
        title('Muscle force');
        saveas(gcf,'results/muscleforce2.fig');
        saveas(gcf,'results/muscleforce2.png');

        plot(abs(buttonheight.data),abs(push_muscle_force.data),'Color','black');
        title('Button displacement vs. Muscle Force');
        saveas(gcf,'results/buttonheight-force.fig');
        saveas(gcf,'results/buttonheight-force.png');

        plot(abs(push_muscle_force.data), abs(buttonheight.data),'Color','black');
        title('Muscle force vs. Button height');
        saveas(gcf,'results/force-buttonheight.fig');
        saveas(gcf,'results/force-buttonheight.png');

        plot(acceleration.time, acceleration.data,'Color','black');
        title('Finger tip acceleration');
        saveas(gcf,'results/tipacceleration.fig');
        saveas(gcf,'results/tipacceleration.png');

        plot(buttonheight.time, buttonheight.data,'Color','black');
        title('Button height');
        saveas(gcf,'results/buttonheight.fig');
        saveas(gcf,'results/buttonheight.png');

        histogram(activation_points);
        title('Make signal points');
        saveas(gcf,'results/activationpoints.fig');
        saveas(gcf,'results/activationpoints.png');

        histogram(percepts);
        title('Integrated p-center');
        saveas(gcf,'results/pcenters.fig');
        saveas(gcf,'results/pcenters.png');

        histogram(peak_ps);
        title('Tactile: p-center');
        saveas(gcf,'results/pcenters-pulp.fig');
        saveas(gcf,'results/pcenters-pulp.png');

        histogram(peak_ds);
        title('Displacement p-center');
        saveas(gcf,'results/pcenters-displacement.fig');
        saveas(gcf,'results/pcenters-displacement.png');

        histogram(peak_ls);
        title('Light: p-center');
        saveas(gcf,'results/pcenters-light.fig');
        saveas(gcf,'results/pcenters-light.png');

        histogram(peak_as);
        title('Audition: p-center');
        saveas(gcf,'results/pcenters-audition.fig');
        saveas(gcf,'results/pcenters-audition.png');

        histogram(muscleforces);
        title('Muscle forces');
        saveas(gcf,'results/muscleforces.fig');
        saveas(gcf,'results/muscleforces.png');

        plot(raw_pulpcontact.time, raw_pulpcontact.data,'Color','black');
        title('Pulp contraction');
        saveas(gcf,'results/pulpcontraction-raw.fig');
        saveas(gcf,'results/pulpcontraction-raw.png');

        plot(raw_displacement.time, raw_displacement.data,'Color','black');
        title('Displacement of finger tip');
        saveas(gcf,'results/displacement-raw.fig');
        saveas(gcf,'results/displacement-raw.png');

        plot(Finger.fstart_height-raw_displacement.data, -velocity.data,'Color','black');
        title('Finger Displacement - Finger Velocity');
        saveas(gcf,'results/displacement-fingervelocity.fig');
        saveas(gcf,'results/displacement-fingervelocity.png');

        plot(displacements, velocities,'Color','black');
        title('Finger Displacement - Finger Velocity');
        saveas(gcf,'results/displacements-fingervelocities.fig');
        saveas(gcf,'results/displacements-fingervelocities.png');

        plot(-buttonheight.data, -velocity.data,'Color','black');
        title('Button Displacement - Finger Velocity');
        saveas(gcf,'results/buttonheight-fingervelocity.fig');
        saveas(gcf,'results/buttonheight-fingervelocity.png');

        plot(-buttonheight.data, raw_pulpcontact.data,'Color','black');
        title('Button Displacement - Pulp Contact');
        saveas(gcf,'results/buttonheight-contact.fig');
        saveas(gcf,'results/buttonheight-contact.png');

    end

% Compute objective function, weighing temporal precision (ovs) and muscle
% effort (muscleforces)
objective = mean(ovs) + mean(muscleforces) * 0.33;
    delta_avg    delta_std    success_avg    ovs_avg    activation_point_avg    activation_point_std    muscleforce_avg    muscleforce_std
    _________    _________    ___________    _______    ____________________    ____________________    _______________    _______________

    0.042828     0.037211     0.89           0.18283    0.37882                 0.053356                1.5215             0.14323        

end
|===============================================================================================================================|
| Iter | Eval   | Objective  | Objective  | BestSoFar  | BestSoFar  |         amag |         pwag |         onag |       exp_ap |
|      | result |            | runtime    | (observed) | (estim.)   |              |              |              |              |
|===============================================================================================================================|
|    1 | Best   |     1.5101 |     28.053 |     1.5101 |     1.5101 |      0.62135 |       3.9435 |       1.4683 |       1.1433 |
|    2 | Best   |     1.2912 |     23.016 |     1.2912 |      1.304 |      0.92765 |       2.0363 |      0.39225 |      0.92023 |
|    3 | Accept |     2.3816 |     29.939 |     1.2912 |     1.2914 |      0.77267 |       4.6403 |       1.0797 |       2.7476 |
|    4 | Accept |     1.3757 |     22.071 |     1.2912 |     1.2912 |      0.80379 |       1.2142 |       1.3966 |       1.8495 |
|    5 | Best   |      1.181 |     25.611 |      1.181 |     1.1812 |        0.524 |       4.6128 |      0.43232 |       0.9853 |
|    6 | Accept |     1.4922 |     29.079 |      1.181 |      1.224 |      0.69704 |       4.3845 |       1.2058 |       2.0619 |
|    7 | Accept |     1.7298 |     23.304 |      1.181 |     1.2741 |       0.6156 |       2.8311 |     0.040666 |       1.4207 |
|    8 | Accept |     2.1406 |     30.983 |      1.181 |     1.1811 |      0.94485 |       4.8265 |       1.4336 |      0.67495 |
|    9 | Best   |     1.0098 |     26.406 |     1.0098 |     1.0099 |       0.5053 |       4.3037 |       1.0636 |      0.88208 |
|   10 | Best   |     1.0055 |     26.935 |     1.0055 |     1.0055 |      0.50136 |       4.4366 |      0.82979 |      0.71387 |
|   11 | Best   |    0.85567 |     23.612 |    0.85567 |    0.85595 |      0.50139 |       3.0941 |      0.74097 |      0.74259 |
|   12 | Accept |     1.3523 |     20.826 |    0.85567 |     0.9276 |      0.50458 |       1.0581 |       0.9277 |       1.2649 |
|   13 | Best   |    0.58859 |     25.244 |    0.58859 |    0.58884 |      0.50513 |       3.7613 |      0.38235 |      0.58384 |
|   14 | Accept |     1.5703 |      25.52 |    0.58859 |    0.61075 |      0.51382 |       3.6367 |      0.89072 |      0.21783 |
|   15 | Accept |    0.67466 |     23.792 |    0.58859 |    0.61312 |      0.50147 |        2.925 |      0.11741 |      0.53586 |
|   16 | Accept |     1.3718 |     20.668 |    0.58859 |    0.58881 |      0.50996 |       1.2075 |      0.29124 |      0.62779 |
|   17 | Accept |     1.6621 |     25.161 |    0.58859 |    0.58875 |      0.55286 |       3.4309 |       1.4429 |      0.57366 |
|   18 | Accept |    0.97351 |     23.076 |    0.58859 |    0.69442 |      0.51071 |       3.7767 |     0.078104 |      0.60784 |
|   19 | Best   |    0.55796 |     24.837 |    0.55796 |    0.55789 |      0.50339 |       3.2367 |      0.39307 |      0.60073 |
|   20 | Accept |    0.60196 |     24.425 |    0.55796 |    0.55788 |       0.5103 |       3.0939 |      0.30016 |      0.45079 |
|===============================================================================================================================|
| Iter | Eval   | Objective  | Objective  | BestSoFar  | BestSoFar  |         amag |         pwag |         onag |       exp_ap |
|      | result |            | runtime    | (observed) | (estim.)   |              |              |              |              |
|===============================================================================================================================|
|   21 | Accept |    0.91476 |     27.464 |    0.55796 |    0.55789 |      0.61524 |        3.336 |      0.36457 |      0.56458 |
|   22 | Accept |    0.62221 |     24.715 |    0.55796 |    0.56668 |      0.50297 |       3.4723 |      0.40586 |        0.535 |
|   23 | Accept |     0.8462 |     23.732 |    0.55796 |    0.60918 |      0.50059 |       3.1965 |      0.38377 |      0.65138 |
|   24 | Accept |    0.59794 |     24.656 |    0.55796 |    0.59779 |      0.50587 |       3.3621 |      0.32111 |      0.47521 |
|   25 | Accept |    0.76113 |     23.868 |    0.55796 |    0.64589 |      0.50943 |       3.3007 |      0.35384 |      0.49728 |
|   26 | Best   |    0.54811 |     25.468 |    0.54811 |    0.61965 |      0.50056 |       3.1227 |      0.18727 |      0.40118 |
|   27 | Accept |    0.72499 |     24.187 |    0.54811 |    0.62565 |      0.50212 |       3.3477 |      0.15103 |      0.33918 |
|   28 | Accept |    0.67325 |     23.019 |    0.54811 |    0.62804 |      0.50122 |       3.1438 |      0.27783 |      0.45813 |
|   29 | Accept |    0.72552 |     23.588 |    0.54811 |    0.63754 |      0.50231 |       3.2461 |      0.27185 |      0.39601 |
|   30 | Accept |     1.1332 |     22.109 |    0.54811 |    0.69021 |      0.50078 |       3.1128 |       0.1986 |      0.43216 |
|   31 | Accept |    0.98932 |     26.083 |    0.54811 |    0.72266 |      0.50008 |       3.4942 |      0.39821 |      0.55219 |
|   32 | Accept |    0.72789 |     23.888 |    0.54811 |    0.72194 |       0.5007 |       3.6619 |      0.30718 |       0.4846 |
|   33 | Accept |    0.81813 |     24.677 |    0.54811 |    0.72792 |      0.50058 |        3.751 |      0.31178 |      0.47442 |
|   34 | Accept |    0.76416 |      25.91 |    0.54811 |    0.72996 |       0.5016 |       3.1983 |      0.34111 |       0.5167 |
|   35 | Accept |    0.57772 |     24.393 |    0.54811 |    0.71755 |      0.50234 |       3.4561 |      0.29292 |      0.52174 |
|   36 | Accept |    0.63157 |     24.598 |    0.54811 |    0.71039 |      0.50464 |       3.8245 |      0.26293 |      0.45875 |
|   37 | Accept |     0.6124 |     26.777 |    0.54811 |    0.70644 |      0.50304 |       3.8594 |      0.18422 |      0.40886 |
|   38 | Accept |    0.65475 |     25.519 |    0.54811 |    0.69091 |      0.50619 |        4.188 |      0.17162 |      0.33661 |
|   39 | Accept |    0.63075 |     26.477 |    0.54811 |    0.68631 |       0.5014 |       4.2682 |      0.17223 |      0.35386 |
|   40 | Accept |    0.65671 |     25.726 |    0.54811 |    0.67262 |      0.50177 |       4.6173 |     0.099236 |      0.27529 |
|===============================================================================================================================|
| Iter | Eval   | Objective  | Objective  | BestSoFar  | BestSoFar  |         amag |         pwag |         onag |       exp_ap |
|      | result |            | runtime    | (observed) | (estim.)   |              |              |              |              |
|===============================================================================================================================|
|   41 | Accept |    0.96019 |      28.17 |    0.54811 |    0.70263 |      0.50039 |       4.6286 |      0.25635 |      0.43321 |
|   42 | Accept |      1.128 |     25.608 |    0.54811 |      0.695 |      0.50006 |       3.8349 |     0.037558 |      0.23365 |
|   43 | Accept |     0.5995 |      27.01 |    0.54811 |    0.68749 |      0.50034 |       3.8529 |      0.31292 |      0.49788 |
|   44 | Accept |    0.59574 |     25.496 |    0.54811 |    0.68129 |      0.50035 |       4.0879 |      0.24637 |      0.47807 |
|   45 | Accept |    0.65959 |     26.799 |    0.54811 |    0.67945 |      0.50255 |       4.0689 |      0.23858 |      0.45165 |
|   46 | Accept |     0.6729 |     27.937 |    0.54811 |    0.67469 |      0.50533 |       4.1576 |       0.2548 |      0.44915 |
|   47 | Accept |    0.65099 |     25.382 |    0.54811 |    0.67219 |      0.50116 |       4.0818 |      0.34767 |      0.50448 |

__________________________________________________________
Optimization completed.
MaxTime of 1200 seconds reached.
Total function evaluations: 47
Total elapsed time: 1215.7282 seconds.
Total objective function evaluation time: 1185.8141

Best observed feasible point:
     amag       pwag      onag      exp_ap 
    _______    ______    _______    _______

    0.50056    3.1227    0.18727    0.40118

Observed objective function value = 0.54811
Estimated objective function value = 0.67219
Function evaluation time = 25.4683

Best estimated feasible point (according to models):
     amag       pwag      onag      exp_ap 
    _______    ______    _______    _______

    0.50464    3.8245    0.26293    0.45875

Estimated objective function value = 0.67219
Estimated function evaluation time = 25.462