You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
How to find the best combinations of calculation of function with 36 input constants, each constant is possible to calculate in 3 ways
5 views (last 30 days)
Show older comments
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
9 Comments
Bruno Luong
on 15 Nov 2023
Edited: Bruno Luong
on 15 Nov 2023
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
5 billions years to evaluate all the combinations.
Still want to do brute force?
Lubos Brezina
on 16 Nov 2023
Ofcourse no. Therefore I asked here. Would you have any advice how to solve it in more efficiently way please?
Bruno Luong
on 16 Nov 2023
"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
on 16 Nov 2023
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
on 17 Nov 2023
Moved: Bruno Luong
on 17 Nov 2023
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
on 17 Nov 2023
If I was you I use ga rather than gamultiobj, by defining the least-squares fit of data to model.
Lubos Brezina
on 18 Nov 2023
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
on 18 Nov 2023
Edited: Bruno Luong
on 18 Nov 2023
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
on 19 Nov 2023
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 :))
Answers (1)
John D'Errico
on 15 Nov 2023
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
4 Comments
Lubos Brezina
on 16 Nov 2023
Edited: Lubos Brezina
on 16 Nov 2023
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
on 16 Nov 2023
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
on 16 Nov 2023
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
on 16 Nov 2023
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.
See Also
Categories
Find more on Surrogate Optimization in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)