Error using fmincon MATLAB function
Show older comments
I would like to obtain the minimum of the function using MATLAB 'fmincon' function where the first eigenvalue derived from the given matrix is my objective function. I was to write/supply its objective function, but it is too long to put it inside the code. I wrote the following code, but I got stuck to continue. The error I found is:
Error using symengine
Unable to convert expression into double array.
I was wondering if I could get help? Thanks.
options = optimoptions('fmincon','GradObj','on');
x = fmincon(@myfun,[1,1],[],[],[],[],[1,1],[2,2],[],options)
function [f,g] = myfun(x)
syms x1 x2;
% 4 by 4 matrix with two variables inside
C=[1 2 x1 4;...
1/2 1 3 x2;...
1/x1 5 1 3/2;...
1/4 1/x2 2 1];
f1 = eig(C); % calculates eigenvalues
f=f1(1); % the first eigenvalue
if nargout > 1 % gradient required
% Calculate partial derivatives of f
g(1) = diff(f,x1);
g(2) = diff(f,x2);
g = double([g(1);g(2)]);% data type double
end
end
2 Comments
Walter Roberson
on 5 Oct 2020
When you say "first" do you mean the one with largest absolute value?
Walter Roberson
on 5 Oct 2020
using the eigenvalue with the largest absolute value is potentially discontinuous in the variables. fmincon requires differentiable jacobian but anything equivalent to a max() operation is not differentiateable at the boundary (and the jacobian potentially turns into zero with respect to some of the variables.)
Accepted Answer
Because of the discontinuity at x1=0 and x2=0, you need to optimize over each quadrant separately.
s=[+1,+1; %quadrant signs
-1,+1;
+1,-1;
-1,-1];
[X1,X2] =ndgrid(linspace(1e-5,2,100));
clear Q xcell
for i=4:-1:1 %optimize over each quadrant
F=arrayfun(@myfun,X1*s(i,1),X2*s(i,2)); %initial discrete search for initial guess
[r,c]=find(F==min(F(:)),1);
x0=[X1(r,c)*s(i,1),X2(r,c)*s(i,2)];
[xcell{i},Q(i)]= fminsearch(@(x)myfun(x(1),x(2)), x0 ,...
optimset('MaxFunEvals',1e8,'MaxIter',1e6));
end
[~,j]=min(Q); %pick the best quadrant
[fOptimal,xOptimal]=deal(Q(j),xcell{j}) %final solution
function fval = myfun(x1,x2)
if abs(abs(x1)-1)>=1, fval=inf; return; end
if abs(abs(x2)-1)>=1, fval=inf; return; end
fval= max(abs(eig([1 2 x1 4;...
1/2 1 3 x2;...
1/x1 5 1 3/2;...
1/4 1/x2 2 1])));
end
8 Comments
No, the solution I gave you is constrained. The bounds given in your post are -2<=x1<=2 and -2<=x2<=2 and that's what the code implements via the lines,
if abs(abs(x1)-1)>=1, fval=inf; return; end
if abs(abs(x2)-1)>=1, fval=inf; return; end
s=[+1,+1]; %quadrant signs
[X1,X2] =ndgrid(linspace(1,2,100));
clear Q xcell
for i=size(s,1):-1:1 %optimize over each quadrant
F=arrayfun(@myfun,X1*s(i,1),X2*s(i,2)); %initial discrete search for initial guess
[r,c]=find(F==min(F(:)),1);
x0=[X1(r,c)*s(i,1),X2(r,c)*s(i,2)];
[xcell{i},Q(i)]= fminsearch(@(x)myfun(x(1),x(2)), x0 ,...
optimset('MaxFunEvals',1e8,'MaxIter',1e6));
end
[~,j]=min(Q); %pick the best quadrant
[fOptimal,xOptimal]=deal(Q(j),xcell{j}) %final solution
function fval = myfun(x1,x2)
if abs(abs(x1)-3/2)>=1/2, fval=inf; return; end
if abs(abs(x2)-3/2)>=1/2, fval=inf; return; end
fval= max(abs(eig([1 2 x1 4;...
1/2 1 3 x2;...
1/x1 5 1 3/2;...
1/4 1/x2 2 1])));
end
Matt J
on 30 Oct 2020
As Walter told you, fmincon is not applicable to this problem because you do not have a differentiable cost function.
Matt J
on 31 Oct 2020
On what intervals are x1,x2 bounded in this case?
HAT
on 31 Oct 2020
Well, I don't know why you think you need to supply your own gradient with fmincon. Using fmincon with its default finite difference approximation produces a minimum in close agreement both with fminbnd and with a surf plot of the objective,
%% Solution via fmincon
[X1,X2]=ndgrid(linspace(1,2,101));
F=arrayfun(@objfunc,X1,X2);
[i,j]=find(F==min(F(:)),1);
x0=[X1(i,j),X2(i,j)];
fun=@(x)objfunc(x(1),x(2));
xfmincon=fmincon(fun,x0, [],[],[],[],[1,1],[+2,+2])
%% Solution via fminbnd
x1=fminbnd(@fun1,1,2);
x2=fminbnd( @(x2)objfunc(x1,x2), 1,2);
xfminbnd=[x1,x2]
%% Visual check
surf(X1,X2,F,'FaceAlpha',.4), xlabel x1, ylabel x2
function f=fun1(x1)
[~,f]=fminbnd( @(x2)objfunc(x1,x2), 1,2);
end
function fval=objfunc(x1,x2)
C=[1 2 exp(x1) 4;...
1/2 1 3 exp(x2);...
1/exp(x1) 1/3 1 1/5;...
4 1/exp(x2) 5 1];
fval=max(abs(eig(C)));
end
HAT
on 1 Nov 2020
Matt J
on 1 Nov 2020
You're welcome, but please Accept-click the answer if you consider the question resolved.
More Answers (1)
Walter Roberson
on 5 Oct 2020
0 votes
You can use the symbolic toolbox to generate eig() on the matrix with Symbolic x1 and x2, and then you can matlabFunction that and pass the handle into the objective function. You would also calculate the g functions with respect to each of the four eigenvalues, and pass those handles as well.
The objective would then not have any symbolic operations. It would call the stored handle upon the input x values, getting back a vector of four eigenvalues. It would figure out which had the largest absolute value and store in f. If nargout requires, it would use the index to select the jacobian handle and evaluate it on the current inputs and return those.
1 Comment
Walter Roberson
on 29 Oct 2020
% fmincon with gradient vector
syms x1 x2
vars = [x1, x2];
C=[1 2 x1 4;...
1/2 1 3 x2;...
1/x1 1/3 1 1/5;...
4 1/x2 5 1];
f=eig(C) % generate eigenvalues on the matrix with Symbolic x1 and x2
eigFunction = matlabFunction(f, 'vars', {vars}); % matlabFunction that and pass the handle into the objective function
% % calculate the g functions with respect to each of the four eigenvalues
g = jacobian([f(1), f(2), f(3),f(4)], vars); % jacobian
gradFunction = matlabFunction(g, 'vars', {vars});
x=ones(2,1); % initial value
% MATLAB fmincon
options = optimoptions('fmincon','GradObj','on');
x = fmincon(@(x)myfun(x, eigFunction, gradFunction), x, [], [], [], [], [1,1], [2,2], [], options)
function [objFunction2 grad2] = myfun(x, eigFunction, gradFunction)
t = eigFunction(x);
[~, maxidx] = max(abs(t));
objFunction2 = t(maxidx);
% If nargout requires, it would use the index
%to select the jacobian handle and evaluate it on the current inputs and return those.
if nargout > 1
gradvalues = gradFunction(x);
grad2 = gradvalues(maxidx);
end
end
Categories
Find more on Calculus in Help Center and File Exchange
