Problem in solving an optimization problem

Hi, guys. This is my code for my optimization problem:

Objective=@Case;
p0 =[1,1,1,1,1];
A = [];
b = [];
Aeq = [];
beq = [];
lb=[1;1;1;1;1];
ub=[7;7;7;7;7];
nonlc=[];
%options = optimoptions('fmincon','Algorithm','interior-point','OptimalityTolerance',1e-20,'FunctionTolerance',1e-20,'StepTolerance',1e-20);

C = fmincon(Objective, p0, A, b, Aeq, beq, lb, ub,nonlc);
disp(C)

function E=Case(p)
moleWt=[28;56;84;112;140;168;156];              % mole weight C2,C4,...,C12,C11 [g/mol]
Woo_Datas_g = {[0.377984 0.696286 2.9244 4.43634 4.63528 4.83422 4.47613],[0.656499 1.69098 5.5504 11.1605 10.882 9.33024 10.6034],[0.69628 1.53183 4.43634 5.74934 6.78382 4.66446 6.98276],[0.497347 1.17374 2.68568 3.6008 2.9244 3.95889 4.07825],[0.497347 1.01459 1.57162 2.2878 2.88462 2.96419 2.48674]};     %Woo's data set for C2,C4,..,C12
F_G_out_mmol = [0.354 0.207 0.169 0.0724 0.0179 0.278 0.291];  %Gas outlet for each case[mmol/s]
t_max = 18000;        %Reaction duration[s]
F_G_out_mol = zeros(1,7); %Gas outlet for each case at the end of the reaction[mol]

for j=1:7
    F_G_out_mol(j) = F_G_out_mmol(j)*t_max*1e-3;
end

C4_mol = zeros(1,7);                              %C4H8 product(mol)
C6_mol = zeros(1,7);                              %C6H12 product(mol)
C8_mol = zeros(1,7);                              %C8H16 product(mol)
C10_mol = zeros(1,7);                             %C10H24 product(mol)
C12_mol = zeros(1,7);                             %C12H224 product(mol)

C4_g = cell2mat(Woo_Datas_g(1));                  %C4H8 product(mol)
C6_g = cell2mat(Woo_Datas_g(2));                  %C6H12 product(mol)
C8_g = cell2mat(Woo_Datas_g(3));                  %C8H16 product(mol)
C10_g = cell2mat(Woo_Datas_g(4));                 %C10H20 product(mol)
C12_g = cell2mat(Woo_Datas_g(5));                 %C12H24 product(mol)

for j=1:7
    C4_mol(j) = C4_g(j)./moleWt(2);
    C6_mol(j) = C6_g(j)./moleWt(3);
    C8_mol(j) = C8_g(j)./moleWt(4);
    C10_mol(j) = C10_g(j)./moleWt(5);
    C12_mol(j) = C12_g(j)./moleWt(6);
end
a = p(1);
b = p(2);
c = p(3);
d = p(4);
e = p(5);
E = (F_G_out_mol(1)-(C4_mol(round(a))+C6_mol(round(b))+C8_mol(round(c))+C10_mol(round(d))+C12_mol(round(e))))^2;
end

As you can see I want E as my objective function. The problem is my parameters for this problem. a,b,c,d,e should be positive integers and they can't be something like 3.5, 0.5 etc. They can only be 1,2,3,..,7. The problem is I can't figure out how define this in fmincon(). Any suggestions for solving this issue?? I have tried to use floor() and run before the parameters but fmincon() gave me a weird answer (1.99 for all the parameters). I also want E to be as minimum as possible. How can I adjust the error in fmincon()??

[Torsten]

I don't understand the sense of your code.

[naiva saeedia]

@Torsten

I completely understand your reason. Sum of the max(Ci_moles) is so small and there is no way it's going to be equal to F_G_out_mol so what's I'm doing is just useless since F_G_out_mol is just greater but again running this code helped me to understand this.

Answer 1

[Manikanta Aditya]

Hi @naiva saeedia

The fmincon function in MATLAB is designed for continuous optimization problems, and it doesn’t directly support integer constraints. However, you can use the ga function from the Global Optimization Toolbox, which supports integer constraints.

Here’s how you can modify your code:

Objective=@Case;
p0 =[1,1,1,1,1,1,1];
A = [];
b = [];
Aeq = [];
beq = [];
lb=[1;1;1;1;1;1;1];
ub=[7;7;7;7;7;7;7];
IntCon = 1:5; % Define the indices of the variables to be integer
options = optimoptions('ga','Display','iter'); % Display iteration information
C = ga(Objective,7,A,b,Aeq,beq,lb,ub,[],IntCon,options);
disp(C)

function E=Case(p)
% ... rest of your code ...
end

In this code, IntCon is a vector that specifies which variables are integer variables, and options is used to set the options for the ga function. The ga function then solves the problem with the integer constraints.

Hope this helps.

[Manikanta Aditya]

Full script can look like this:

Objective=@Case;
p0 =[1,1,1,1,1,1,1];
A = [];
b = [];
Aeq = [];
beq = [];
lb=[1;1;1;1;1;1;1];
ub=[7;7;7;7;7;7;7];
IntCon = 1:5; % Define the indices of the variables to be integer
options = optimoptions('ga','Display','iter'); % Display iteration information
C = ga(Objective,7,A,b,Aeq,beq,lb,ub,[],IntCon,options);
disp(C)

function E = Case(p)
    moleWt = [28; 56; 84; 112; 140; 168; 156]; % mole weight C2,C4,...,C12,C11 [g/mol]
    Woo_Datas_g = {[0.377984 0.696286 2.9244 4.43634 4.63528 4.83422 4.47613], [0.656499 1.69098 5.5504 11.1605 10.882 9.33024 10.6034], [0.69628 1.53183 4.43634 5.74934 6.78382 4.66446 6.98276], [0.497347 1.17374 2.68568 3.6008 2.9244 3.95889 4.07825], [0.497347 1.01459 1.57162 2.2878 2.88462 2.96419 2.48674]}; % Woo's data set for C2,C4,..,C12
    F_G_out_mmol = [0.354 0.207 0.169 0.0724 0.0179 0.278 0.291]; % Gas outlet for each case [mmol/s]
    t_max = 18000; % Reaction duration [s]
    F_G_out_mol = zeros(1, 7); % Gas outlet for each case at the end of the reaction [mol]

    for j = 1:7
        F_G_out_mol(j) = F_G_out_mmol(j) * t_max * 1e-3;
    end

    C4_mol = zeros(1, 7); % C4H8 product (mol)
    C6_mol = zeros(1, 7); % C6H12 product (mol)
    C8_mol = zeros(1, 7); % C8H16 product (mol)
    C10_mol = zeros(1, 7); % C10H24 product (mol)
    C12_mol = zeros(1, 7); % C12H224 product (mol)

    C4_g = cell2mat(Woo_Datas_g(1)); % C4H8 product (mol)
    C6_g = cell2mat(Woo_Datas_g(2)); % C6H12 product (mol)
    C8_g = cell2mat(Woo_Datas_g(3)); % C8H16 product (mol)
    C10_g = cell2mat(Woo_Datas_g(4)); % C10H20 product (mol)
    C12_g = cell2mat(Woo_Datas_g(5)); % C12H24 product (mol)

    for j = 1:7
        C4_mol(j) = C4_g(j) ./ moleWt(2);
        C6_mol(j) = C6_g(j) ./ moleWt(3);
        C8_mol(j) = C8_g(j) ./ moleWt(4);
        C10_mol(j) = C10_g(j) ./ moleWt(5);
        C12_mol(j) = C12_g(j) ./ moleWt(6);
    end

    a = round(p(1));
    b = round(p(2));
    c = round(p(3));
    d = round(p(4));
    e = round(p(5));

    E = (F_G_out_mol(1) - (C4_mol(a) + C6_mol(b) + C8_mol(c) + C10_mol(d) + C12_mol(e)))^2;
end

[Sam Chak]

Hi @naiva saeedia

In the static optimization problem, it is assumed that the data is measured with ultra-precision exactly at 18000 seconds. Although a minimum can be determined, it is important to evaluate the significance of the cost function value "fval = 36.5262". Additionally, you need to establish the maximum value of fval that you would consider reliable. A value of "fval = 0" indicates that the perfect solution has been found.

Objective = @Case;
p0        = [1,1,1,1,1,1,1];
A         = [];
b         = [];
Aeq       = [];
beq       = [];
lb        = [1;1;1;1;1;1;1];
ub        = [7;7;7;7;7;7;7];
IntCon    = 1:5; % Define the indices of the variables to be integer
[C, fval] = ga(Objective,7,A,b,Aeq,beq,lb,ub,[],IntCon)

ga stopped because the average change in the penalty function value is less than options.FunctionTolerance and 
the constraint violation is less than options.ConstraintTolerance.

C = 1x7

6.0000 4.0000 7.0000 7.0000 6.0000 3.1897 3.6924

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fval = 36.5262

function E = Case(p)
    moleWt      = [28; 56; 84; 112; 140; 168; 156]; % mole weight C2,C4,...,C12,C11 [g/mol]
    Woo_Datas_g = {[0.377984 0.696286 2.9244 4.43634 4.63528 4.83422 4.47613], [0.656499 1.69098 5.5504 11.1605 10.882 9.33024 10.6034], [0.69628 1.53183 4.43634 5.74934 6.78382 4.66446 6.98276], [0.497347 1.17374 2.68568 3.6008 2.9244 3.95889 4.07825], [0.497347 1.01459 1.57162 2.2878 2.88462 2.96419 2.48674]}; % Woo's data set for C2,C4,..,C12
    F_G_out_mmol= [0.354 0.207 0.169 0.0724 0.0179 0.278 0.291]; % Gas outlet for each case [mmol/s]
    t_max       = 18000; % Reaction duration [s]
    F_G_out_mol = zeros(1, 7); % Gas outlet for each case at the end of the reaction [mol]

    for j = 1:7
        F_G_out_mol(j) = F_G_out_mmol(j) * t_max * 1e-3;
    end

    C4_mol  = zeros(1, 7); % C4H8 product (mol)
    C6_mol  = zeros(1, 7); % C6H12 product (mol)
    C8_mol  = zeros(1, 7); % C8H16 product (mol)
    C10_mol = zeros(1, 7); % C10H24 product (mol)
    C12_mol = zeros(1, 7); % C12H224 product (mol)

    C4_g    = cell2mat(Woo_Datas_g(1)); % C4H8 product (mol)
    C6_g    = cell2mat(Woo_Datas_g(2)); % C6H12 product (mol)
    C8_g    = cell2mat(Woo_Datas_g(3)); % C8H16 product (mol)
    C10_g   = cell2mat(Woo_Datas_g(4)); % C10H20 product (mol)
    C12_g   = cell2mat(Woo_Datas_g(5)); % C12H24 product (mol)

    for j = 1:7
        C4_mol(j)   = C4_g(j)  ./ moleWt(2);
        C6_mol(j)   = C6_g(j)  ./ moleWt(3);
        C8_mol(j)   = C8_g(j)  ./ moleWt(4);
        C10_mol(j)  = C10_g(j) ./ moleWt(5);
        C12_mol(j)  = C12_g(j) ./ moleWt(6);
    end

    a = round(p(1));
    b = round(p(2));
    c = round(p(3));
    d = round(p(4));
    e = round(p(5));

    E = (F_G_out_mol(1) - (C4_mol(a) + C6_mol(b) + C8_mol(c) + C10_mol(d) + C12_mol(e)))^2;
end

[naiva saeedia]

@Manikanta Aditya

Thanks a lot for your help. It has solved my problem.

[naiva saeedia]

@Sam Chak

Yeah exactly for the perfect solution I need fval to be zero but since my F_G_out_mol(1) is greater than sum of the maximum of(Ci_mol) It's never happening in my code. I guess this just shows F_G_out is not equal to sum of the Ci_mol.

[Manikanta Aditya]

Thanks @naiva saeedia

[naiva saeedia]

@Manikanta Aditya

I have a question about ga() optimizer. I have seven cases in my problems as you can see(Seven arrays in F_G_out_mol Vector) I want to run the optimization for each case.(In my code here I only run it for the first case F_G_out_mol(1)) Matlab gave me these answers for example:

a =1,b=2,c=2,d=4,e=6

When I run the optimization problem for my second case(Changing F_G_out_mol(1) to F_G_out_mol(2)) I want solver to choose a value for a from 2,3,4,5,6,7, choose a value for b from 1,3,4,5,6,7, choose value for c from 1,3,4,5,6,7, choose value for d from 1,2,3,5,6,7, etc. I guess you understood my pattern here. I want to delete the privious answer that ga() gave me an find unique values for a,b,c,d,e for each case. Is this possible?

[Manikanta Aditya]

Hi @naiva saeedia, To solve the optimization problem for each case with the requirement that the genetic algorithm (ga) chooses values for (a, b, c, d, e) without repeating the values used in the previous case, you can modify the bounds or constraints dynamically based on the results of the previous optimization.

However, MATLAB's ga function does not natively support excluding specific values directly within its standard options. Instead, you can achieve this by customizing the fitness function or by adjusting the bounds for each run based on the results from the previous optimization.

Objective = @Case;
A = [];
b = [];
Aeq = [];
beq = [];
IntCon = 1:5; % Define the indices of the variables to be integer

% Initialize variables to store the results
results = zeros(7, 5); % Assuming 7 cases and 5 decision variables

for i = 1:7
    if i == 1
        % For the first case, use the initial bounds
        lb = [1;1;1;1;1];
        ub = [7;7;7;7;7];
    else
        % For subsequent cases, adjust bounds based on previous results
        lb = ones(5, 1);
        ub = 7 * ones(5, 1);
        
        % Note: This is a conceptual approach; you'll need to adjust it based on your actual constraints
        for j = 1:length(IntCon)
            prevValue = results(i-1, j);
            % Here, you might need a more sophisticated way to adjust bounds or add custom constraints
            % For example, you might adjust `lb` and `ub` to exclude `prevValue`
        end
    end
    
    % Update the Objective function to pass the current case index
    currentObjective = @(p) Case(p, i);
    
    % Run the genetic algorithm
    options = optimoptions('ga', 'Display', 'iter', 'InitialPopulationMatrix', lb'); % Assuming you adjust options as needed
    [C, fval] = ga(currentObjective, 5, A, b, Aeq, beq, lb, ub, [], IntCon, options);
    
    % Store the results
    results(i, :) = C;
    
    % Display or process results as needed
    disp(['Case ', num2str(i), ':']);
    disp(C);
end

function E = Case(p, caseIndex)
    % Your fitness function here, modified to use `caseIndex` if necessary
end

Hope this clarifies.

[Sam Chak]

Hi Naiva, could you please formulate the problem using easy-to-understand mathematical equations (not in code form)? This will enable Mani to guide you more effectively.

For instance, if there are equations with unknown coefficients, it is believed that there are certain values of the coefficients that would result in a balance between the equations and the given data.

[naiva saeedia]

@Manikanta Aditya

Thanks a lot for sharing your concept. That's what I have done by hand(Trying the first case then adjust the boundaries for the next case and so on). I have tried to adjust your code but I think what I want is not that simple. For example when I first ran your code without adjusting boundaries I have got these results:

Case1:2,3,3,3,2

Case2:1,1,6,4,6

Case3:1,2,2,7,6

so on. As you can see case3 have repeated answers.

I have tried to change the case3 lb and ub in your for loop. I have kept the boundary conditions for case 1 and 2 the same since they don't produce the repeated answers like case3. I have tried these values for lb and up:

lb=[3;4;4;5;3]

ub=[7;7;7;7;7]

now as you can see we still have problems for lb and ub in case3. For example for the second array(ub and lb of b parameter) the 2 is exculded from ub and lb but 2 should be there! and only 1 and 3 shouldn't be in ub and lb. for third array(ub and lb of c parameter) 6 still in the ub and lb and 1 and 2 exculded. I really don't know how should I solve this problem. I guess you have told me to use prevValue to solve this problem. but how can I add and remove specific values from ub and lb in each case. I don't understand this part. I guess we just can't do this. Am I wrong??

[naiva saeedia]

Hi @Sam Chak. The concept of the code is simple I guess. I have five matricies:

C4_mol

C6_mol

C8_mol

C10_mol

C12_mol

These are my products in the reactor. I also have one matrix that shows the products that produced (F_G_out). We ran our Reactor 7 times. Each time is a case test. So I have seven arrays for each of these matricies. I don't know which array in Ci_mol matricies belongs to which case. If we use mass balance sum of the products should be equal to F_G_out. I want to find out which array in Ci_mol belongs to which case and that's why I have defined this optimization problem. I want to find the indicies of C4_mol,C6_mol,..,C12_mol matricies that gave me the F_G_out(1). It's the same for other cases.

[Torsten]

Make 5 nested loops from 1 to 7, call your function with the respective 5 parameters and compute E. This will give you 7^5 results for E. Order them in ascending order and extract the one for which indices don't repeat.

[naiva saeedia]

@Torsten

I guess I didn't quite understand what you are trying to say. In this way I get results for some cases and for some parameters. I want to find the results without indices repeatation for every case and every parameter.

[Torsten]

Since you can only work with integer values from 1 to 7 for the 5 parameters (since they are array indices), the 7^5 cases are all cases you need to consider, aren't they ?

[naiva saeedia]

@Torsten

I need to calculate all of the parameters for seven cases that means I need to calculate 35 variables in my code. If what you meant is what I have desceribed then yes(I didn't quite understand what you meant by 7^5, why you have considered so many cases??)

[Torsten]

p(1) can have values 1,2,3,4,5,6,7

p(2) can have values 1,2,3,4,5,6,7

...

p(5) can have values 1,2,3,4,5,6,7

->

[p(1),p(2),..,p(5)] can be 7*7*7*7*7 = 7^5 different vectors.

At least this is what your code suggests.

[naiva saeedia]

@Torsten

Sorry for late reply. Yes you are correct. Thanks for explanation about 7^5 vectors.