How to find the best combinations of calculation of function with 36 input constants, each constant is possible to calculate in 3 ways
2 次查看(过去 30 天)
显示 更早的评论
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 个评论
Bruno Luong
2023-11-15
编辑:Bruno Luong
2023-11-15
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
2023-11-16
Ofcourse no. Therefore I asked here. Would you have any advice how to solve it in more efficiently way please?
Bruno Luong
2023-11-16
"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
2023-11-16
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
2023-11-17
移动:Bruno Luong
2023-11-17
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
2023-11-17
If I was you I use ga rather than gamultiobj, by defining the least-squares fit of data to model.
Lubos Brezina
2023-11-18
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
2023-11-18
编辑:Bruno Luong
2023-11-18
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
2023-11-19
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 :))
回答(1 个)
John D'Errico
2023-11-15
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 个评论
Lubos Brezina
2023-11-16
编辑:Lubos Brezina
2023-11-16
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
2023-11-16
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
2023-11-16
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
2023-11-16
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Parallel Computing Fundamentals 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- 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)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)