Why is my gamultobj algorithm not working?

Dear MATLAB community,

I have a question regarding a MOO algorithm I am currently implementing in order to optimize a charging current profile for my thesis. For the optimization I am using gamultobj. My fitness function is quite complex as you can see in the code example below, it includes calculations for three different ode's and a time loop over 900 timesteps.

My problem occurs when I try to add a Nonlinear fitness function into my problem using nonlcon. The algorithm just doesn't finish calculating when I include the nonlcon function into my algorithm (even after long waiting time). Without it, it runs through in acceptable time interval. Also, it doesn't throw out any error messages, the algorithm just continues on calculating.

My nonlinear constraint function basically does the same calculations as my fitness function, it just calculates the constraints at the end instead of the fitness values.

Is there anything I am missing about the nonlinear constraint functions or do you have any other clue what I can do to debug and optimize my code? I would be grateful for any help!

I attached a simplified version of my code including the fitness and nonlinear function below.

% parameters
population = 2;
generations = 2;
convergence = 1e-3;
paretoFraction = 0.35;

options = optimoptions('gamultiobj', ...
    'PopulationSize', population, ...
    'MaxGenerations', generations, ...
    'FunctionTolerance', convergence, ...
    'ParetoFraction', paretoFraction, ...
    'UseParallel',true,...
    'UseVectorized', false,...
    'OutputFcn', @nsga2_outputfcn, ...
    'PlotFcn', @gaplotpareto);

fitnessFunction = @(x) multiObjectiveFunction_Multistage(x);
nonlcon = @(x) constraintFunction(x);
[x, fval, exitflag, output, population, scores] = gamultiobj(fitnessFunction, nVars, [], [], [], [], lb, ub, nonlcon, options);

function objectives = multiObjectiveFunction_Multistage(x)
    global start_SOC Q_cell_As nStages_max stagelength COP_Heizen COP_Kuehlen T_initial cell_voltage_initial delta_SOC ntimesteps

    % weighting factors
    omega = [0.5, 0.1, 0.4];

    % initialisation
    f2 = zeros(1, ntimesteps);
    f3 = zeros(1, ntimesteps);

    % overpotentials
    U_RC_1_c_i = 0;
    U_RC_2_c_i = 0;
    U_RC_1_a_i = 0;
    U_RC_2_a_i = 0;

    % cell voltage
    U_cell_last = cell_voltage_initial;

    I_ch = zeros(1, nStages_max);
    I_heat = zeros(1, nStages_max);

    for i = 1:nStages_max
        I_ch(i) = x(2 * i - 1);
        I_heat(i) = x(2* i);
    end

    T = T_initial;

    stage = 1;
    SOC_x = start_SOC;

    %% loop

    while stage <= nStages_max && SOC_x < start_SOC + delta_SOC

        I_ch_stage = I_ch(stage);
        I_heat_stage = I_heat(stage);

        t = 1;

        [t1a, U_RC1_a_vec, t2a, U_RC2_a_vec] = overpotential_a_v3(I_ch_stage,SOC_x,U_RC_1_c_i,U_RC_2_c_i,T);
        [t1c, U_RC1_c_vec, t2c, U_RC2_c_vec] = overpotential_c_v3(I_ch_stage,SOC_x,U_RC_1_c_i,U_RC_2_c_i,T);

        [tT, T_vec] = zelltemperatur_v3(I_ch_stage, I_heat_stage, U_RC_1_a_i, U_RC_2_a_i, U_RC_1_c_i, U_RC_2_c_i, T, U_cell_last);

        while t <= stagelength && SOC_x < start_SOC + delta_SOC

            if stage == 1 && t == 1
                SOC_x = start_SOC;
            else
                SOC_x = SOC_x + (I_ch_stage / Q_cell_As);
            end

            U_RC1_a = interp1(t1a, U_RC1_a_vec, t);
            U_RC2_a = interp1(t2a, U_RC2_a_vec, t);
            U_RC1_c = interp1(t1c, U_RC1_c_vec, t);
            U_RC2_c = interp1(t2c, U_RC2_c_vec, t);

            T = interp1(tT, T_vec, t);

            [parameters_c, parameters_a] = interpolation_v3(SOC_x,T);
            U_ocv_a = parameters_a(1);
            U_ocv_c = parameters_c(1);
            interp_R_0_a = parameters_a(2);
            interp_R_0_c = parameters_c(2);

            % overpotenials
            eta_a = I_ch_stage * interp_R_0_a + U_RC1_a + U_RC2_a; 
            eta_c = I_ch_stage * interp_R_0_c + U_RC1_c + U_RC2_c;

            % ocv
            U_ocv_cell = U_ocv_c - U_ocv_a;

            % electrode voltages
            U_c = U_ocv_c + eta_c;
            U_a = U_ocv_a - eta_a;

            % cell voltage
            U_cell = U_c - U_a;

            % calculation f2
            if I_heat_stage < 0
                f2((stage - 1) * stagelength + t) = I_ch_stage * (U_cell - U_ocv_cell) + I_heat_stage * U_cell * COP_Kuehlen;
            elseif I_heat_stage >= 0
                f2((stage - 1) * stagelength + t) = I_ch_stage * (U_cell - U_ocv_cell) + I_heat_stage * U_cell * COP_Heizen;
            end

            % calculation f3
            f3((stage - 1) * stagelength + t) = platingpower_v3(I_ch_stage, SOC_x, U_RC1_a, U_RC2_a, T, interp_R_0_a);

            % Update cell voltage
            U_cell_last = U_cell;

            % timestep
            t = t + 1;

        end

        U_RC_1_c_i = U_RC1_c;
        U_RC_2_c_i = U_RC2_c;
        U_RC_1_a_i = U_RC1_a;
        U_RC_2_a_i = U_RC2_a;

        stage = stage + 1;

    end

    total_f1 = omega(1) * (delta_SOC * Q_cell_As)/mean(I_ch);
    total_f2 = omega(2) * sum(f2);
    total_f3 = omega(3) * sum(f3);

    objectives = [total_f1, total_f2, total_f3];
end

function [c, ceq] = constraintFunction_vereinfacht_2(x)

    global ntimesteps nStages_max stagelength T_initial start_SOC delta_SOC Q_cell_As cell_voltage_initial T_max voltage_max voltage_min T_min
    
    U_RC_1_c_i = 0;
    U_RC_2_c_i = 0;
    U_RC_1_a_i = 0;
    U_RC_2_a_i = 0;
    
    c = zeros(ntimesteps,4);
    
    T = T_initial;
    U_cell_last = cell_voltage_initial;
    SOC_x = start_SOC;
    
    stage = 1;

    while stage <= nStages_max && SOC_x < start_SOC + delta_SOC

        I_ch_stage = x(2 * stage - 1);
        I_heat_stage = x(2 * stage);

        t = 1;

        [t1a, U_RC1_a_vec, t2a, U_RC2_a_vec] = overpotential_a_v3(I_ch_stage,SOC_x,U_RC_1_c_i,U_RC_2_c_i,T);
        [t1c, U_RC1_c_vec, t2c, U_RC2_c_vec] = overpotential_c_v3(I_ch_stage,SOC_x,U_RC_1_c_i,U_RC_2_c_i,T);

        [tT, T_vec] = zelltemperatur_v3(I_ch_stage, I_heat_stage, U_RC_1_a_i, U_RC_2_a_i, U_RC_1_c_i, U_RC_2_c_i, T, U_cell_last);

        while t <= stagelength && SOC_x < start_SOC + delta_SOC
            
            if stage == 1 && t == 1
                SOC_x = start_SOC;
            else
                SOC_x = SOC_x + (I_ch_stage / Q_cell_As);
            end

            U_RC1_a = interp1(t1a, U_RC1_a_vec, t);
            U_RC2_a = interp1(t2a, U_RC2_a_vec, t);
            U_RC1_c = interp1(t1c, U_RC1_c_vec, t);
            U_RC2_c = interp1(t2c, U_RC2_c_vec, t);

            T = interp1(tT, T_vec, t);
        
            [parameters_c, parameters_a] = interpolation_v3(SOC_x,T);
            U_ocv_a = parameters_a(1);
            U_ocv_c = parameters_c(1);
            interp_R_0_a = parameters_a(2);
            interp_R_0_c = parameters_c(2);

            eta_a = I_ch_stage * interp_R_0_a + U_RC1_a + U_RC2_a; 
            eta_c = I_ch_stage * interp_R_0_c + U_RC1_c + U_RC2_c;

            U_c = U_ocv_c + eta_c;
            U_a = U_ocv_a - eta_a;

            U_cell = U_c - U_a;

            % Constraint 1: maximum temperature
            c((stage - 1) * stagelength + t,1) =  T - T_max;

            % Constraint 2: minimum temperature
            c((stage - 1) * stagelength + t,2) =  5 - T_min;

            % Constraint 3: maximum cell voltage
            c((stage - 1) * stagelength + t,3) =  U_cell - voltage_max;

            % Constraint 4: minimum cell voltage
            c((stage - 1) * stagelength + t,4) =  voltage_min - U_cell;

            U_cell_last = U_cell;

        end

        U_RC_1_c_i = U_RC1_c;
        U_RC_2_c_i = U_RC2_c;
        U_RC_1_a_i = U_RC1_a;
        U_RC_2_a_i = U_RC2_a;

        stage = stage + 1;

    end
 
    ceq = [];

end

Answer 1

[Matt J]

My guess would be that it is struggling to find feasible solutions to the nonlinear constraints (perhaps because they don't have any). One way to investigate this would be to run the code with the objective function replaced by a trivial single objective,

fitnessFunction=@(x) 0;

so that the only thing it needs to do to complete the optimization is to find a feasible x.

[Walter Roberson]

In situations like that, I just start by assuming that one of the global variables was not defined, so one of the expressions involving the variable came out empty and that is why you encounter the problem.

This is not always the solution, but it is the solution more often than not. Don't use globals.

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html