How to find the best combinations of calculation of function with 36 input constants, each constant is possible to calculate in 3 ways

Hi community,

I have a function with 36 constants (it possible to image the function as an surface). Each of the 36 constant is possible to calculate in 3 ways (named/substituted for example as 0, 1, 2). It is totally 3 ^ 36 of possible combinations. I am looking for the best combination to find the minimum sum error from reference values set.

My idea to solve it was with brute force at first: for each combination 1 to 3 ^ 36 to calculate input control vector in ternary system (0...000, 0...001, 0...002, 0...010, etc.) and then to evaluate function output (based on control vector) for each reference values set. If the solution was better than previous one, I stored the "better combination" and go to the next combination. But the solution takes unreal time.

It is there any better solution how to solve it? I read something about optimalization toolbox, but I am not familiar with it. I read also about parfor function, but my solution described above is evaluating every step if is better than previous one. So, I don’t know, if it is even possible to calculate it in parallel way.

Thank you in Advance

[Bruno Luong]

Let assume you can evalute 1000 combinations by second it will take you

nopspersec = 1000;
nbsecperh = 3600;
nbhperday = 24;
nbdayperyear = 365;
timeinyear = 3^36/nopspersec/nbsecperh/nbhperday/nbdayperyear

timeinyear = 4.7595e+06

5 billions years to evaluate all the combinations.

Still want to do brute force?

[Lubos Brezina]

Ofcourse no. Therefore I asked here. Would you have any advice how to solve it in more efficiently way please?

[Bruno Luong]

"Would you have any advice how to solve it in more efficiently way please?"

You only describe your problem very surperficially (36 variables each take 3 value). The dependency of the objectiive to those parameters can be then anything. There is no way we can give you any useful advide from that little information.

[Steven Lord]

As an aside, an easier way to compute how long an operation will take (at N combinations per second) is to use duration functionality.

format longg
numCombinations = 3^36

numCombinations = 
      1.50094635296999e+17

combinationsPerSecond = 1000;
yearsRequired = years(seconds(numCombinations/combinationsPerSecond))

yearsRequired = 
          4756309.64920183

[Lubos Brezina]

I appologise. I needed more time to prepare some aditional information and to study GA, how it works. Basically, target is to fit orange outputs of my function to blue reference points with minimum error (deviation). It is like surface fitting by control of 36 integer values, which control a mixing of 36 constants (named as "Ainputs").

I try to us gamultiobj: each reference point is as separate object. Main body is:

%%Preparation and imports
clear;

%%Import of A inputs (36x2)
opts = spreadsheetImportOptions("NumVariables", 2);
File =pwd +"\Ainput.xlsx";
opts.Sheet = "List1";
opts.VariableTypes = ["double","double"];
opts.DataRange = "A1:B36";
Ainput = table2array(readtable(File, opts, "UseExcel", false));
clear opts File

%%Import of reference values x_ref, y_ref, value_ref
File =pwd +"\RefValues.xlsx";
opts = spreadsheetImportOptions("NumVariables", 1);
opts.Sheet = "List1";
opts.VariableTypes = "double";
opts.MissingRule = "omitrow";
% Import the data
opts.DataRange = "A2:A10000";
x_ref= table2array(readtable(File, opts, "UseExcel", false));
opts.DataRange = "B2:B10000";
y_ref = table2array(readtable(File, opts, "UseExcel", false));
opts.DataRange = "C2:C10000";
value_ref = table2array(readtable(File, opts, "UseExcel", false));
clear opts File

NoOfRefPoints = numel(value_ref);
MixingRatio = [0.5, 0.5];
NoOfUnits = 2;

%%Default values before optimalization - just for visualization of problem
defMixingControl = [0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];
DefaultValues = PointDevFunction(defMixingControl, NoOfUnits, MixingRatio, Ainput,NoOfRefPoints, x_ref, y_ref, value_ref);
hScatter = scatter3([x_ref, x_ref],[y_ref, y_ref],[value_ref, DefaultValues']);
legend(hScatter, 'Value_ref', 'Default values', 'Location', 'NorthEast', 'Interpreter', 'none' );
xlabel('x');
ylabel('y');
zlabel('z');

%%gamultiobj preparation and run
nVars=36;
intcon=1:nVars;
rng default
A=[];
b=[];
Aeq=[];
beq=[];
%MixingControl based on 0 or 1 or 2 integer values
lb=[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];
ub=[2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2;2];
nonlcon=[];  

fun = @(MixingControl)PointDevFunction(MixingControl, NoOfUnits, MixingRatio, Ainput,NoOfRefPoints, x_ref, y_ref, value_ref);
optsM = optimoptions(@gamultiobj, 'PlotFcn',{@gaplotdistance,@gaplotpareto,@gaplotrankhist, @gaplotstopping,@gaplotspread});    
optsM.MaxGenerations=200*36;
optsM.MaxStallGenerations=1E5;
SavedOptions = optsM;
[x,FvalmexitFlag,Output]=gamultiobj(fun,nVars,A,b,Aeq,beq,lb,ub,nonlcon,intcon, optsM);

PointDevFunction do firstly once mixing of "A" constants based on MixingControl matrix (36x1), where each control cell is integer 0 or 1 or 2. As next step is done main calculation for every reference point (like surface) and obtain a value and its PointDeviation: which should be as minimum as possible.

function PointDeviation = PointDevFunction(MixingControl, NoOfUnits, Ratio, Ainput, NoOfRefPoints, x, y, value_ref)
    
% Mixing calculation    
Amixed = MixingFunction(NoOfUnits, Ratio, Ainput, MixingControl); 
    
% Calcualtion of points deviation in terms of reference points x(ii), y(ii), value_ref
PointDeviation(1:NoOfRefPoints)=0;
for ii = 1:NoOfRefPoints
    %main function
    B2 = Amixed(1) * x(ii)^2 + Amixed(2) * x(ii) + Amixed(3) + Amixed(4) / x(ii) + Amixed(5) / x(ii)^2 + Amixed(6) / x(ii)^3 + Amixed(7) / x(ii)^4 + Amixed(8) / x(ii)^5 + Amixed(9) / x(ii)^6;
    B3 = Amixed(10) * x(ii)^2 + Amixed(11) * x(ii) + Amixed(12) + Amixed(13) / x(ii) + Amixed(14) / x(ii)^2 + Amixed(15) / x(ii)^3 + Amixed(16) / x(ii)^4 + Amixed(17) / x(ii)^5 + Amixed(18) / x(ii)^6;
    B4 = Amixed(19) + Amixed(20) / x(ii) + Amixed(21) / x(ii)^2;
    B5 = Amixed(22) + Amixed(23) / x(ii) + Amixed(24) / x(ii)^2;
    B6 = Amixed(25) + Amixed(26) / x(ii) + Amixed(27) / x(ii)^2;
    B7 = Amixed(28) + Amixed(29) / x(ii) + Amixed(30) / x(ii)^2;
    B8 = Amixed(31) + Amixed(32) / x(ii) + Amixed(33) / x(ii)^2;
    B9 = Amixed(34) + Amixed(35) * x(ii) + Amixed(36) / x(ii);
    value = 1 + B2 * y(ii) + B3 * y(ii)^2 + B4 * y(ii)^3 + B5 * y(ii)^4 + B6 * y(ii)^5 + y(ii)^2 / x(ii)^3 * (B7 + B8 * y(ii)^2) * (exp(-B9 * y(ii)^2));    

    %deviation
    PointDeviation(ii)=value - value_ref(ii);
end        

MixingFunction combine 3 differents mixing rules for every "A" constants (36x2), according to MixingControl vector (36x1):

function Amixed = MixingFunction(NoOfUnits, ratio, Ainput, MixingControl)
    % mixing function, based on MixingControl matrix (0, 1, 2)
    Amixed(1:36) = 0;
    for jj = 1:36        
        if MixingControl(jj) == 0
            for ii = 1:NoOfUnits
                Amixed(jj) = Amixed(jj) + ratio(ii) * Ainput(jj, ii); %suma x*I
            end            
        elseif MixingControl(jj) == 1
            for ii = 1:NoOfUnits                
                Amixed(jj) = Amixed(jj) + ratio(ii) * Ainput(jj, ii) / sqrt(abs(Ainput(jj, ii))); %            
            end
            Amixed(jj) = Amixed(jj) * abs(Amixed(jj));            
        elseif MixingControl(jj) == 2
            for ii = 1:NoOfUnits
                Amixed(jj) = Amixed(jj) + ratio(ii) * Ainput(jj, ii) / ((Ainput(jj, ii) ^ 2) ^ (1/3));
            end
            Amixed(jj) = Amixed(jj) ^ 3;
        else
            error('Non defined MixingControl value')
        end
    end
end

Yes, my syntax is not a proffesional one...I am still learning.

I don´t know, if I am using gamultiobj in good way. It is like stucked now without any progress. Output x (mixing control vector) is also 4x36. Why? I would expect 1x36 only.

Thank you, LB

[Bruno Luong]

If I was you I use ga rather than gamultiobj, by defining the least-squares fit of data to model.

[Lubos Brezina]

Thanks Bruno. I used it at first (ga). Output value for ga was sum of errors of all reference points. Is this +- it what you mean? Ga was without any progress after 4 hours run. Maybe the least square could help. How would you do implementation of "by defining the least-squares fit of data to model"? I have no idea now (I always use it only in passive form - f.e. in curve fitter toolbox).

[Bruno Luong]

A least squares problem is minimizing the objective function defined as

norm(abs(model(param) - data)^2

or

sum(abs(model(param) - data).^2)

sum is over output variables (your surface data).

Such objective function usually eases the minimizer for convergence.

Might be you problem is simply too hard to solve

[Lubos Brezina]

A least square implementation significantly increase quality of results (ga). It was good advice. Ga founded some solutions. Question is, if there could be any better solutions. Because, it was quite hard to force ga to not be stalled at beginning. After hours of "debugging", I changed Crossovefraction to 0.6 (I tried change anything) and it helped. (solution after 30 sec). But that was all. This solution is not the best, I expected more. Does it have any meaning to keep ga running, f.e. 1 day more? I dont think so.

Gamultiobj is stalled. I let it go.

With brute force, I validated ga result. It was possible, because solution was at "side" of 36 constants range (only 18 positions of 36). Other "side" of control matrix (the rest) would take 37years of calculation :))

Answer 1

[John D'Errico]

Brute force is NEVER going to be a good idea. Instead, if you have the global optimization toolbox...

help ga

 GA    Constrained optimization using genetic algorithm.
    GA attempts to solve problems of the following forms:
        min F(X)  subject to:  A*X  <= B, Aeq*X  = Beq (linear constraints)
         X                     C(X) <= 0, Ceq(X) = 0 (nonlinear constraints)
                               LB <= X <= UB 
                               X(i) integer, where i is in the index
                               vector INTCON (integer constraints)
 
    Note: If INTCON is not empty, then no equality constraints are allowed.
    That is:-
    * Aeq and Beq must be empty
    * Ceq returned from NONLCON must be empty
 
    X = GA(FITNESSFCN,NVARS) finds a local unconstrained minimum X to the
    FITNESSFCN using GA. NVARS is the dimension (number of design
    variables) of the FITNESSFCN. FITNESSFCN accepts a vector X of size
    1-by-NVARS, and returns a scalar evaluated at X.
 
    X = GA(FITNESSFCN,NVARS,A,b) finds a local minimum X to the function
    FITNESSFCN, subject to the linear inequalities A*X <= B. Linear
    constraints are not satisfied when the PopulationType option is set to
    'bitString' or 'custom'. See the documentation for details.
 
    X = GA(FITNESSFCN,NVARS,A,b,Aeq,beq) finds a local minimum X to the
    function FITNESSFCN, subject to the linear equalities Aeq*X = beq as
    well as A*X <= B. (Set A=[] and B=[] if no inequalities exist.) Linear
    constraints are not satisfied when the PopulationType option is set to
    'bitString' or 'custom'. See the documentation for details.
 
    X = GA(FITNESSFCN,NVARS,A,b,Aeq,beq,lb,ub) defines a set of lower and
    upper bounds on the design variables, X, so that a solution is found in
    the range lb <= X <= ub. Use empty matrices for lb and ub if no bounds
    exist. Set lb(i) = -Inf if X(i) is unbounded below;  set ub(i) = Inf if
    X(i) is unbounded above. Linear constraints are not satisfied when the
    PopulationType option is set to 'bitString' or 'custom'. See the 
    documentation for details.
 
    X = GA(FITNESSFCN,NVARS,A,b,Aeq,beq,lb,ub,NONLCON) subjects the
    minimization to the constraints defined in NONLCON. The function
    NONLCON accepts X and returns the vectors C and Ceq, representing the
    nonlinear inequalities and equalities respectively. GA minimizes
    FITNESSFCN such that C(X)<=0 and Ceq(X)=0. (Set lb=[] and/or ub=[] if
    no bounds exist.) Nonlinear constraints are not satisfied when the
    PopulationType option is set to 'bitString' or 'custom'. See the 
    documentation for details.
 
    X = GA(FITNESSFCN,NVARS,A,b,Aeq,beq,lb,ub,NONLCON,options) minimizes
    with the default optimization parameters replaced by values in OPTIONS.
    OPTIONS can be created with the OPTIMOPTIONS function. See OPTIMOPTIONS
    for details. For a list of options accepted by GA refer to the
    documentation.
 
    X = GA(FITNESSFCN,NVARS,A,b,[],[],lb,ub,NONLCON,INTCON) requires that
    the variables listed in INTCON take integer values. Note that GA does
    not solve problems with integer and equality constraints. Pass empty
    matrices for the Aeq and beq inputs if INTCON is not empty.
 
    X = GA(FITNESSFCN,NVARS,A,b,[],[],lb,ub,NONLCON,INTCON,options)
    minimizes with integer constraints and the default optimization
    parameters replaced by values in OPTIONS. OPTIONS can be created with
    the OPTIMOPTIONS function. See OPTIMOPTIONS for details.
 
    X = GA(PROBLEM) finds the minimum for PROBLEM. PROBLEM is a structure
    that has the following fields:
        fitnessfcn: <Fitness function>
             nvars: <Number of design variables>
             Aineq: <A matrix for inequality constraints>
             bineq: <b vector for inequality constraints>
               Aeq: <Aeq matrix for equality constraints>
               beq: <beq vector for equality constraints>
                lb: <Lower bound on X>
                ub: <Upper bound on X>
           nonlcon: <Nonlinear constraint function>
            intcon: <Index vector for integer variables>
           options: <Options created with optimoptions('ga',...)>
          rngstate: <State of the random number generator>
 
    [X,FVAL] = GA(FITNESSFCN, ...) returns FVAL, the value of the fitness
    function FITNESSFCN at the solution X.
 
    [X,FVAL,EXITFLAG] = GA(FITNESSFCN, ...) returns EXITFLAG which
    describes the exit condition of GA. Possible values of EXITFLAG and the
    corresponding exit conditions are
 
      1 Average change in value of the fitness function over
         options.MaxStallGenerations generations less than 
         options.FunctionTolerance and constraint violation less than 
         options.ConstraintTolerance.
      3 The value of the fitness function did not change in
         options.MaxStallGenerations generations and constraint violation 
         less than options.ConstraintTolerance.
      4 Magnitude of step smaller than machine precision and constraint
         violation less than options.ConstraintTolerance. This exit 
         condition applies only to nonlinear constraints.
      5 Fitness limit reached and constraint violation less than
         options.ConstraintTolerance. 
      0 Maximum number of generations exceeded.
     -1 Optimization terminated by the output or plot function.
     -2 No feasible point found.
     -4 Stall time limit exceeded.
     -5 Time limit exceeded.
 
    [X,FVAL,EXITFLAG,OUTPUT] = GA(FITNESSFCN, ...) returns a
    structure OUTPUT with the following information:
              rngstate: <State of the random number generator before GA started>
           generations: <Total generations, excluding HybridFcn iterations>
             funccount: <Total function evaluations>
         maxconstraint: <Maximum constraint violation>, if any
               message: <GA termination message>
 
    [X,FVAL,EXITFLAG,OUTPUT,POPULATION] = GA(FITNESSFCN, ...) returns the
    final POPULATION at termination.
 
    [X,FVAL,EXITFLAG,OUTPUT,POPULATION,SCORES] = GA(FITNESSFCN, ...) returns
    the SCORES of the final POPULATION.
 
 
    Example:
      Unconstrained minimization of Rastrigins function:
        function scores = myRastriginsFcn(pop)
            scores = 10.0 * size(pop,2) + sum(pop.^2 - 10.0*cos(2*pi .* pop),2);
 
      numberOfVariables = 2
      x = ga(@myRastriginsFcn,numberOfVariables)
 
      Display plotting functions while GA minimizes
       options = optimoptions('ga','PlotFcn',...
         {@gaplotbestf,@gaplotbestindiv,@gaplotexpectation,@gaplotstopping});
       [x,fval,exitflag,output] = ga(fitfcn,2,[],[],[],[],[],[],[],options)
 
    An example with inequality constraints and lower bounds
     A = [1 1; -1 2; 2 1];  b = [2; 2; 3];  lb = zeros(2,1);
     fitfcn = @(x)0.5*x(1)^2 + x(2)^2 -x(1)*x(2) -2*x(1) - 6.0*x(2);
 
     % Use mutation function which can handle constraints
     options = optimoptions('ga','MutationFcn',@mutationadaptfeasible);
     [x,fval,exitflag] = ga(fitfcn,2,A,b,[],[],lb,[],[],options);
 
    If FITNESSFCN or NONLCON are parameterized, you can use anonymous
    functions to capture the problem-dependent parameters. Suppose you want
    to minimize the fitness given in the function myfit, subject to the
    nonlinear constraint myconstr, where these two functions are
    parameterized by their second argument a1 and a2, respectively. Here
    myfit and myconstr are MATLAB file functions such as
 
         function f = myfit(x,a1)
         f = exp(x(1))*(4*x(1)^2 + 2*x(2)^2 + 4*x(1)*x(2) + 2*x(2) + a1);
 
    and
 
         function [c,ceq] = myconstr(x,a2)
         c = [1.5 + x(1)*x(2) - x(1) - x(2);
               -x(1)*x(2) - a2];
         % No nonlinear equality constraints:
          ceq = [];
 
    To optimize for specific values of a1 and a2, first assign the values
    to these two parameters. Then create two one-argument anonymous
    functions that capture the values of a1 and a2, and call myfit and
    myconstr with two arguments. Finally, pass these anonymous functions to
    GA:
 
      a1 = 1; a2 = 10; % define parameters first
      % Mutation function for constrained minimization
      options = optimoptions('ga','MutationFcn',@mutationadaptfeasible);
      x = ga(@(x)myfit(x,a1),2,[],[],[],[],[],[],@(x)myconstr(x,a2),options)
 
    Example: Solving a mixed-integer optimization problem
    An example of optimizing a function where a subset of the variables are
    required to be integers:
    
    % Define the objective and call GA. Here variables x(2) and x(3) will
    % be integer. 
    fun = @(x) (x(1) - 0.2)^2 + (x(2) - 1.7)^2 + (x(3) -5.1)^2; 
    x = ga(fun,3,[],[],[],[],[],[],[],[2 3])
           
    See also OPTIMOPTIONS, FITNESSFUNCTION, GAOUTPUTFCNTEMPLATE, PATTERNSEARCH, @.

    Documentation for ga
       doc ga

[Lubos Brezina]

Thank you. I try it (i have never used GA before). If you can advice me on which parts of GA I need to focus regarding to my problem, It would be appreciated.

EDIT: deleted

[John D'Errico]

Which parts of GA? You just need to read the help for GA. Look at the examples. I'm not sure what else to say.

You have 36 variables, all of which will take on 3 different levels, so you can think of them as having the levels 1,2, or 3. Then use those as an index to generate the actual levels of those variables.

[Torsten]

@John D'Errico

I have no knowledge about "ga". Maybe you can help here: Is "ga" really better suited for this problem than simply testing a fixed number of random triples ?

[Lubos Brezina]

Basically: I want to fit surface given by my function to reference surface (reference points) by mentioned 36 constants (1, 2, 3).

Based on your recommendation, I tried it with GA. My attitude is, that function output is sum of errors of all reference points. It works, but I don´t know, if it converts/leads to any solution.

I am also trying GAmultiobj. Each error of each reference point is a separate output of function. It doesn´t seems to be work well now, but maybe it need to be tuned more.

P.S. f.e. these I meant about focusing: If it is better for GA or GAmultiobj, etc.