# Find gradient of objective function for fmincon

5 views (last 30 days)

Show older comments

Arnab Samaddar-Chaudhuri
on 11 Jan 2023

Commented: Arnab Samaddar-Chaudhuri
on 12 Jan 2023

Hello,

I need to optimize 2 points. I am finding it difficult to manually calculate the gradient of the objective function The shortened version of the code is as follows:

W1 = @(X) objfun(X,XQ,YQ,cdf1,cdf2);

% I have a meshgrid, where XQ, YQ, cdf1, cdf2 are of mxn matrix and are known

options = optimoptions('fmincon','Algorithm','sqp', 'Display','iter-detailed',...

'CheckGradients',true,'SpecifyObjectiveGradient',true);

lb = [1909,-390;

2059,-390];

ub = [1969,-360;

2059,-360];

initValue = [1969,-360;

2059,-360];

%these values are not important here

[x1,~,exitflag] = fmincon(W1,initValue,[],[],[],[], lb1, ub1, [], options);

%where

function [f,gradf] = objfun(X,XQ,YQ,cdf1,cdf2)

p1 = @(Point) interp2(XQ,YQ,cdf1,Point(1,1),Point(1,2),'nearest');

p2 = @(Point) interp2(XQ,YQ,cdf2,Point(2,1),Point(2,2),'nearest');

f = -1+(1-p1(X))*(1-p2(X));

if nargout > 1

[gx1, gy1] = gradient(cdf1);

[gx2, gy2] = gradient(cdf2);

dp1x = @(Point) interp2(XQ,YQ,gx1,Point(1,1),Point(1,2),'nearest');

dp1y = @(Point) interp2(XQ,YQ,gy1,Point(1,1),Point(1,2),'nearest');

dp2x = @(Point) interp2(XQ,YQ,gx2,Point(2,1),Point(2,2),'nearest');

dp2y = @(Point) interp2(XQ,YQ,gy2,Point(2,1),Point(2,2),'nearest');

gradfx1 = @(P) - dp1x * (1-p2(P)) ;

gradfy1 = @(P) - dp1y * (1-p2(P)) ;

gradfx2 = @(P) - (1-p1(P)) * dp2x ;

gradfy2 = @(P) - (1-p1(P)) * dp2y ;

gradfx = [gradfx1 ; gradfx2];

gradfy = [gradfy1 ; gradfy2];

gradf = [gradfx;gradfy];

end

end

I get the following error:

Error using vertcat

Nonscalar arrays of function handles are not allowed; use cell arrays instead.

Error in SBR_Optimisation_tool_1_8_4/objfun (line 2893)

gradfx = [gradfx1 ; gradfx2];

I am sure I am not calculating the gradient correctly. I simpy cannot use gradient(f) or gradient(p1). That gives other kinds of error. How to calculate the gradient here?

##### 4 Comments

Torsten
on 11 Jan 2023

Edited: Torsten
on 11 Jan 2023

p1 = @(x,y) interp2(XQ,YQ,cdf1,x,y,'nearest');

p2 = @(x,y) interp2(XQ,YQ,cdf2,x,y,'nearest');

These two functions are not even continuous functions of x and y. So how could you calculate a gradient for them ?

Further, you must supply numerical values for the gradient at x, not a function handle as you do.

Be happy if fmincon works with your non-differentiable function and don't try to make things even worse by supplying non-existing gradients.

Walter Roberson
on 11 Jan 2023

To get around that particular problem change

gradfx = [gradfx1 ; gradfx2];

gradfy = [gradfy1 ; gradfy2];

gradf = [gradfx;gradfy];

to

gradfx = {gradfx1 ; gradfx2};

gradfy = {gradfy1 ; gradfy2};

gradf = [gradfx;gradfy];

The result will be a 4 x 1 cell array of function handles, which is legal in MATLAB.

MATLAB does not permit you to use [] to put more than one function handle into an array, because if it did, gradfx(1) would be ambiguous about whether you are invoking the array of function handles each on the input [1], or if you were invoking each of the function handles with no input and then indexing the result at the first location, or if you were selecting the first function handle to be returned as a function handle, or if you were selecting the first function handle and invoking it with no input. It therefore forces you to store multiple function handles in cell, as you can then use {} and () to control your meaning.

### Accepted Answer

Matt J
on 11 Jan 2023

Edited: Matt J
on 11 Jan 2023

Perhaps as follows. Note my restructuring of your array into an Nx2x2 form.

function [f,gradf] = objfun(X,XQ,YQ,cdf1,cdf2)

[x1,y1,x2,y2]=deal(X(:,1,1), X(:,2,1), X(:,1,2), X(:,2,2));

p1 = interp2(XQ,YQ,cdf1,x1,y1,'cubic');

p2 = interp2(XQ,YQ,cdf2,x2,y2,'cubic');

f = -1+(1-p1)*(1-p2);

if nargout > 1

D=1e7*eps(X);

[dx1,dy1,dx2,dy2]=deal(D(:,1,1), D(:,2,1), D(:,1,2), D(:,2,2));

p1x1=(interp2(XQ,YQ,cdf1,x1+dx1,y1,'cubic') -p1)./dx1;

p1y1=(interp2(XQ,YQ,cdf1,dx1,y1+dy1,'cubic') -p1)./dy1;

p2x2=(interp2(XQ,YQ,cdf2,x2+dx2,y2,'cubic') -p2)./dx2;

p2y2=(interp2(XQ,YQ,cdf2,dx2,y2+dy2,'cubic') -p2)./dy2;

gradp1=[p1x1;

p1y1].*(1-p2);

gradp2=[p2x2;

p2y2].*(1-p1);

gradf=[gradp1;gradp2];

end

end

##### 2 Comments

### More Answers (1)

Walter Roberson
on 11 Jan 2023

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!