How do I set up fmincon correctly for parallel computation? "Supplied objective function must return a scalar value"
显示 更早的评论
I have an expensive objective function which takes in multiple parameters, while the contraints are constant. It looks something like this:
out1, out2 = obj_func(in1, in2)
in1 has 1000 different cases. For each possible value of in1, I want to minimize out1 by adjusting in2. At the same time, out2 has to be < 1 in all cases. Because obj_func takes a while to run, I want to parallelize the computation, either by calculating several gradients at the same time, or solve the problem for several in1.
Below is my simplified code:
nPoints = length(INVARS(:,1));
odInputNames = {'DISA'; 'ALT'; 'MACH'};
% Number of batches depending on given batch size:
sizeBatch = 20;
nBatch = floor( nPoints/sizeBatch ) + 1; % Last batch contains remaining points
% Define constraints
CONSTR = struct(...
'param_A', '...',...
'param_B', '...',...
'param_C', '...',...
CONSTR_REQ = [ 1800; ...
1.10; ...
1.15; ...
];
for i = 1:nBatch
% Create batchwise input vectors
for iODVar = 1:length(odInputNames)
if i < nBatch
OD.(odInputNames{iODVar}) = INVARS(sizeBatch*(i-1)+1:sizeBatch*i, iODVar);
else
OD.(odInputNames{iODVar}) = INVARS(sizeBatch*(i-1)+1:end, iODVar); % Last batch has fewer points
end
end
% Setup of FMINCON Optimizer
% ---------------------------------------------------------------------
A = [];
b = [];
Aeq = [];
beq = [];
lb = ones(size(OD.DISA)) * 80;
ub = ones(size(OD.DISA)) * 120;
options = optimoptions(@fmincon, ...
'Algorithm', 'active-set',...
'MaxIter', 200,...
'ConstraintTolerance', 1e-8,...
'OptimalityTolerance', 1e-8,...
'Display', 'off',...
'UseParallel', true);
x0 = ones(size(OD.DISA)) * 85;
% Declare constraint results to be shared between obj_fun and con_fun
constr_fields = fieldnames(CONSTR);
con_res = zeros(length(constr_fields), length(OD.DISA));
% Execute Optimizer
[PLA_max, obj_min,exitflag,~,~] = fmincon( @(x)obj_fun(x), x0, A, b, Aeq, beq, lb, ub, @(x)con_fun(), options );
% Objective:
function [obj] = obj_fun( PLA )
OD.PLA = PLA;
[OD.OUT, OD.IN] = expensive_function( DES, OD );
obj = - OD.OUT.param_obj;
for i=1:length(constr_fields)
con_res(i,:) = eval(['OD.' CONSTR.(constr_fields{i})]);
end
end
% Constraints:
function [c,ceq] = con_fun()
exp_CONSTR_REQ = repmat(CONSTR_REQ, 1, length(OD.DISA));
c = (con_res - exp_CONSTR_REQ); % Non-Equality constraint ( !<= 0 )
ceq = [];
end
It works fine when I solve one case at a time, but when I try to parallelize it I got this error. How can I fix it?
Error using fmincon (line 348)
Supplied objective function must return a scalar value.
2 个评论
Why is lb > ub ?
lb = ones(size(OD.DISA)) * 120;
ub = ones(size(OD.DISA)) * 80;
Why does the error message report something about fminunc ? It is not referenced in your code.
Shao-Yu
2025-1-22
采纳的回答
I want to parallelize the computation, either by calculating several gradients at the same time, or solve the problem for several in1.
The UseParallel option, which you are already invoking, does parallelize the gradient computation. If there is an analytical algorithm for computing the gradient, however, that will often be quicker.
The error message about fmincon returning a scalar should be easy for you to investigate just by calling the objective function on the initial point x0. If obj_fun(x0) gives you a non-scalar value, then you have a bug, and you must trace it.
If you are deliberately having it return a non-scalar, then you have a misconception about how fmincon works. fmincon cannot perform multiple optimizations in a single call.
25 个评论
Shao-Yu
2025-1-23
Matt J
2025-1-23
Pose the optimization problem as the sum of the objective functions that you want to mimize
Torsten
2025-1-23
Why do you return
obj = - sum(OD.OUT.param_obj);
and not
obj = sum(OD.OUT.param_obj);
?
Shao-Yu
2025-1-23
OD.PLA = PLA; % PLA is an array
Shouldn't DES be an array as well? In my suggestion, I indicated that there should be a different DES for every term in the objective function,
Are you saying you can't reproduce your old results even when all the DES are set equal to each other des_i = des_j ?
I could only get it to work when all arrays here were single doubles instead.
I still don't know what this means. Originally your x1 was an array of size(OD.DISA). Unless you are saying that OD.DISA was a scalar, then PLA was never a "single double".
Also, it is not clear what was working before compared to now. With your original objective function f(x), what values of in1,in2 gave you a reasonable optimization result? And what happens with your new objective function F(x1,x2,x3,x4) when you set all four pairs of in1,in2 to the same pair that worked for the optimization of f()?
And have you verified that F(x1,x2,x3,x4) returns the same value as f(x1) + f(x2) + f(x3) + f(x4) up to floating point round-offs?That is a basic validation step that should be done before any optimization is attempted.
Torsten
2025-1-24
What happens if you use
in1 = [0 0 0];
in2 = [10 10 10];
?
Shao-Yu
2025-1-24
What happens if you use
function [obj] = obj_fun( PLA )
in1 = OD.in1;
in2 = OD.in2;
tmp = zeros(size(PLA));
for i = 1:numel(PLA)
OD.PLA = PLA(i);
OD.in1 = in1(i);
OD.in2 = in2(i);
[OD.OUT, OD.IN] = expensive_function(OD);
tmp(i) = OD.OUT.param_obj;
end
obj = - sum(tmp);
OD.in1 = in1;
OD.in2 = in2;
% OD.OUT.param_obj is an array of doubles that each correspond to each
% combination of [in1(1), in2(1)], [in1(2), in2(2)] and so on
end
Shao-Yu
2025-1-24
It's just a test whether your "expensive_function" works correctly.
If the code
function [obj] = obj_fun( PLA )
in1 = OD.in1;
in2 = OD.in2;
tmp = zeros(size(PLA));
for i = 1:numel(PLA)
OD.PLA = PLA(i);
OD.in1 = in1(i);
OD.in2 = in2(i);
[OD.OUT, OD.IN] = expensive_function(OD);
tmp(i) = OD.OUT.param_obj;
end
obj = - sum(tmp);
OD.in1 = in1;
OD.in2 = in2;
% OD.OUT.param_obj is an array of doubles that each correspond to each
% combination of [in1(1), in2(1)], [in1(2), in2(2)] and so on
end
works and the code
function [obj] = obj_fun( PLA )
OD.PLA = PLA; % PLA is an array
[OD.OUT, OD.IN] = expensive_function(OD);
obj = - sum(OD.OUT.param_obj);
% OD.OUT.param_obj is an array of doubles that each correspond to each
% combination of [in1(1), in2(1)], [in1(2), in2(2)] and so on
end
does not work, the problem must lie in your "expensive_function".
It seems that evaluation of array inputs to "expensive_function" for in1, in2 and PLA give different results in param_obj as if "expensive_function" were called for each array element separately.
I would test this outside the optimization environment.
Matt J
2025-1-24
I'm glad you got it working, but if you have now achieved a solution through this answer, please Accept-click it.
I only see the option to call fmincon 20 times in a parfor loop for 1 PLA at a time. This doesn't use the option to evaluate your expensive_function for several PLA inputs, but it should be faster than using 20 PLAs in 1 call to "fmincon". But maybe 20 is still too little to see an advantage of the parallel option taking into account the overhead of parallelization.
Ok, it worked but sadly not faster at all, probably because you can't optimize 20 independent variables with only one objective...
Certainly you can. Even without UseParallel, it is unclear why you don't see speed-up based on what you've told us. You yourself claimed that expensive_function() could process 20 PLA's faster as a batch than with 20 separate calls. That alone should have offered some improvement.
Maybe calling "fmincon" 20 times with one solution variable is faster than calling it once with 20 solution variables.
Like solving 20 (independent) ODEs in one call to an ode integrator will most probably be slower than solving 1 ODE at a time in 20 calls to the ode integrator.
Come to think of it, Torsten is right, but primarily because you are not using the SpecifyObjectiveGradient and SpecifyObjectiveConstraints options. If you have a way of supplying the objective and constraint gradients analytically, you should do so, and it should make things faster.
If you have to rely on the default finite difference approach, then batch computations will not help you. fmincon calls the objective 2*N times where N is the number of independent variables. So, if you try to optimize in batches of N PLAs, the number of function calls required per iteration also grows by a factor of N.
If you have a way of analytically supplying the objective and constraint gradients analytically, you should do so, and it should make things faster.
Even if you don't have an analytical formula for the derivatives, you could try to approximate them with your own finite difference computation. In this special case where your objective function is a separable sum, the finite difference derivatives can be computed as
where E() is the vector-valued output of expensive_function(). This results in only two calls to expensive_function() per iteration instead of N. Similar things apply to your nonlinear cnostraint functions.
Shao-Yu
2025-1-24
Seems like the only options left are making better initial guesses for PLA and approximating the derivatives.
And also Torsten's suggestion of using a parfor loop. That might be the simplest.
So you're saying I should make another function to calculate E(x_i) and E(x_i + δ)
The gradient computations must be done inside your objective functions as described here,
Obviously, you are free to out-source parts of that computation with calls to other functions, if you wish.
Should I just set δ to be something very small, or is there a better way to do it?
You will need to experiment with δ. The CheckGradients fmincon option can help you see if the derivatives you've provided agree with fmincon's own numerical approximation.
Shao-Yu
2025-1-27
更多回答（1 个）
Catalytic
2025-1-24
0 个投票
There are other acceleration strategies that may be superior to parallelization.
I notice in your examples that the in1,in2 pairs change by small amounts. You also appear to believe that the final solutions for PLA will be close to each other, otherwise why initialize the optimization with the same x0? Those being the case, it may be worthwile to go back to the strategy of optimizing one PLA at a time. However, you would use the optimized PLA from the N-th optimization as the x0 for the next optimization. This might cut down on the number of iterations you need to get the solution, and reduce the number of calls to expensive function.
类别
在 帮助中心 和 File Exchange 中查找有关 Get Started with Optimization Toolbox 的更多信息
Translated by 
