How do I input linear constraints for the genetic algorithm function using discrete variables?

Hello,

I am using the matlab's ga algorithm to solve an optimization problem using discrete variables. My decision variable is set to take only integer values, then using a simple mapping function, my decision variable are transformed to the n'th indexed value in a list of specific values which i provide, where n is the integer value initially assigned to the decision variable. I followed the exact steps of the following link:

https://fr.mathworks.com/help/gads/solving-a-mixed-integer-engineering-design-problem-using-the-genetic-algorithm.html

Within this link, they inputed non-linear constraints using a function. In my case, I am interested in inserting linear constraints of inequality, using the matrices A and b. The algorithm seems to not respect these constraints. This is possible because the algorithm is applying the constraints to the decision variables in the form of integers rather than their transformed form.

How should I inser the A and b matrices in such a scenario?

Thank you.

[Torsten]

Where do you use the integer decision variables in their transformed form ?

[Samer Alsamadi]

The objective function uses the decision variables in their transformed form. It starts off by calling a function called EVMapVariables(x), then uses the transformed form of the decision variables to calculate a cost.

[Torsten]

If you don't use nonlinear equality constraints together with the integer decision variables in function "nonlcon", I don't know why the procedure you describe should not work.

But of course we cannot give further advice without code that reproduces the behaviour you describe. What MATLAB version are you using ?

[Samer Alsamadi]

I just added the files to this comment. You can run the Main_Code.m which will run the ga algorithm for 10 generations then stop. The linear inequality constraints I provided using matrix A and b are to ensure that the transformed form of the solution is greater or equal to a certain value in the vector t_ar. So solution x(1) >= t_ar(1), x(2) >= t_ar(2) and so on. To do this I coded A as a diagonal matrix of -1's and b as the negative of t_ar, since the constraint assumes the form: A * X <= b.

However, executing the code "xbestDisc >= t_ar" at the end indicates that sometimes the constraint is respected and sometimes its not. Should I somehow present the A and b matrices in another manner? Thanks again.

[Torsten]

Use lb and ub instead of A and b. Then I think the constraints should always be respected during the solution process.

[Samer Alsamadi]

The example I gave above includes only 1 of the various linear inequality constraints I intend to add to the program. Therefore it will not be possible to just use lb and ub. For example 1 of my constraints is: x(2) >= x(1) + c, where c is a constant value. I don't recall reading anywhere that the ga function with discrete variables can not be given linear constraints. Maybe I'm missing something here?

[Torsten]

I don't recall reading anywhere that the ga function with discrete variables can not be given linear constraints. Maybe I'm missing something here?

You can use linear constraints, but I think if you put them in A and b, it's not guaranteed that they are respected during the solution process. The final result should of course respect them, though.

This is different if you use lower and upper bounds: as far as I can remember, they are respected throughout the complete solution process. And I thought this is what you were asking for, weren't you ?

[Samer Alsamadi]

For clarification, the lb and ub limits I put were intended for the integer form of the decision variables, and they were indeed respected.

The A and b matrix constraints however were intended for the transformed form of the decision variables. These were not respected.

So I guess the question here is: how do I tell the algorithm to apply the A and b linear inequality constraints to the transformed version of the decision variables rather than the integer one.

[Walter Roberson]

I don't think you can do what you want in your version, except by enforcing all of the integer constraints yourself after telling ga() that you are not using any integer constraints.

[Torsten]

The A and b matrix constraints however were intended for the transformed form of the decision variables. These were not respected.

It is not possible to constrain the transformed variables with A and b. With A and b, you can only constrain solution variables. But the solution variables are the decision integers, not the deduced transformed variables.

You have to define the nonlinear inequality constraints on the transformed variables in "nonlcon".

But I don't understand why you wrote the constraints look like

x(1) >= t_ar(1), x(2) >= t_ar(2) 

Aren't x(1), x(2) the integer decision variables and t_ar(1), t_ar(2) some fixed constants known right from the beginning of the computation ?

[Torsten]

I just see in your code that you work with random inputs that change during each iteration. This is stochastic optimization and not possible with the usual MATLAB optimizers.

clear all
clc
rng('default') % For reproducibility

N = 10; % Nbr of electric vehicles
eta = 0.9;
C = 50; %kWh
P_ev = 8; %kWh
SoC_i = [0.5, 0.4, 0.6, 0.8, 0.1, 0.2, 0.5, 0.85, 0.55, 0.35];
SoC_f = 1;
time_step = 0.25; % 15 min. planning interval

delta_xi = zeros(1,10);
for i = 1:N
    delta_xi_temp = ((SoC_f - SoC_i(i)) * C) / (eta * P_ev);
    time_step = 0.25;
    delta_xi(i) = ceil(delta_xi_temp / time_step) * time_step;
end

t_ar = [4,5,6,7,8,9,10,11,12,13];
t_dp = t_ar + 8; 


A = diag(-1*ones(1,10)); % Arrival time inequality left side
   
b = [-t_ar(1); 
      -t_ar(2);
      -t_ar(3);
      -t_ar(4);
      -t_ar(5);
      -t_ar(6);
      -t_ar(7);
      -t_ar(8);
      -t_ar(9);
      -t_ar(10)]; % Arrival time inequality right side
    
lb = ones(1,10);
ub = 96*ones(1,10);

nvars = 10;

opts = optimoptions(@ga, ...
                    'PopulationSize', 150, ...
                    'MaxGenerations', 10, ...
                    'EliteCount', 10, ...
                    'FunctionTolerance', 1e-8, ...
                    'ConstraintTolerance',1.0000e-06, ...
                    'PlotFcn', @gaplotbestf);
                
  
[xbestDisc, fbestDisc, exitflag] = ga(@EVWithDisc, ...
    10, A, b, [], [], lb, ub, [], 1:10, opts);

Optimization terminated: maximum number of generations exceeded.

xbestDisc = EVMapVariables(xbestDisc);

xbestDisc >= t_ar

ans = 1×10 logical array
   1   1   0   1   0   0   0   0   0   0

function cost = EVWithDisc(x)
%cantileverVolumeWithDisc Calculate volume of a stepped cantilever
%
%   V = cantileverVolumeWithDisc(x) calculates the volume of the stepped
%   cantilever in the "Solving a Mixed Integer Engineering Design Problem
%   Using the Genetic Algorithm" example. In this function, the variables
%   x(3), x(4), x(5) and x(6) are chosen from a discrete set. 

%   Copyright 2012 The MathWorks, Inc.

% Map the discrete variables
x = EVMapVariables(x);

% Volume of cantilever beam
cost = EV_fitness_func_(x);
end

function x = EVMapVariables(x)

all_X = 0:0.25:23.75;

x(1) = all_X(x(1));
x(2) = all_X(x(2));
x(3) = all_X(x(3));
x(4) = all_X(x(4));
x(5) = all_X(x(5));
x(6) = all_X(x(6));
x(7) = all_X(x(7));
x(8) = all_X(x(8));
x(9) = all_X(x(9));
x(10) = all_X(x(10));

end

function cost = EV_fitness_func_(chargingTimes)
    rng('default')
    a_t = 10 * rand(1,96);
    N = 10; % Nbr of electric vehicles
    eta = 0.9;
    C = 50; %kWh
    P_ev = 8; %kWh
    SoC_i = [0.5, 0.4, 0.6, 0.8, 0.1, 0.2, 0.5, 0.85, 0.55, 0.35];
    SoC_f = 1;
    time_step = 0.25; % 15 minute time steps for planning
    
    cost = 0;
    delta_xi = zeros(1,10);
    for i = 1:N
       x_i = chargingTimes(i);
       sum_internal = 0;
       delta_xi(i) = ((SoC_f - SoC_i(i)) * C) / (eta * P_ev);
       delta_xi(i) = ceil(delta_xi(i) / time_step) * time_step;
       for t = x_i:0.25:x_i+delta_xi(i)
           t_temp = t;
           if t_temp > 23.75
               t_temp = t_temp - 24;
           end
           index = (t_temp/0.25) + 1;
           electricity_price = a_t(index);
           sum_internal = sum_internal + P_ev * electricity_price * 0.25;
       end
       cost = cost + sum_internal;
    end
end

[Samer Alsamadi]

The algorithm's inputs are actually deterministic. The only "random" part of it is a_t, however I did not intend for it to be random. I'm just using " a_t = 10 * rand(1,96) " to generate the matrix a_t.

Recalling that my goal was to implement linear constraints of inequality on the transformed form of my decision variables, I managed to do that by simply transforming my constraints to be implemented on the "integer" form of my variables. The reason for which I called it simple is because the transformation function used to transform the variable from the "integer" form to the "transformed" form is simple. It's a function that replaces the decision variable with the i'th indexed value of a discrete list of values.

Example:

Integer form: x = 5.

Discrete_list = [0,0.25,0.5,0.75,1,1.25,1.5,1.75,2].

Transformed form: x = Discrete_list(5) = 1.

Thank you everyone for your assistance!

[Torsten]

a_t and thus the electricity price change with each call to your cost function. So how should "ga" be able to get a stable deterministic solution to your problem ?

[Samer Alsamadi]

I believe the rng('default') line of code at the start should lead to the "pseudo-random" matrix a_t being the same every time. Am I incorrect in thinking so ?

[Torsten]

Only the seed of the random generator is set to a specified value that remains constant each time you run the code. But if you call the random number generator several times within one program execution, the output changes:

rng("default")
rand()

ans = 0.8147

rand()

ans = 0.9058

rand()

ans = 0.1270

[Samer Alsamadi]

Thanks @Torsten for your clarification, i'll make sure I just specify a_t myself rather than using a random number generator.

[Torsten]

Inequality constraints on the transformed variables have to be set like this:

[xbestDisc, fbestDisc, exitflag] = ga(@EVWithDisc, ...
                                      10, [],[], [], [], lb, ub, @nonlcon, 1:10, opts);

function [c,ceq] = nonlcon(x)
  all_X = 0:0.25:23.75;
  t_ar = [4,5,6,7,8,9,10,11,12,13];
  n = numel(x);
  x_trans(1:n) = all_X(x(1:n));
  c = -x_trans + t_ar;
  ceq = [];
end

[Samer Alsamadi]

I wanted to try that method but read the writing linear constraints in the non-linear constraints part of the algorithm could lead to problems as the algo will mess things up. Do you have an idea if that is true or not? Thanks anyways for providing a code example!

[Torsten]

Your constraints are linear in the transformed variables, but not in the solution variables.

Or are you able to use A (which is multiplied by x, not x_trans) and b to set the constraints you want on x_trans ?

What you did in your code is setting x(1:10) >= t_ar(1:10), not x_trans(1:10) >= t_ar(1:10). And from what you wrote this is what you actually want, isn't it ?