I am looking to find the minimum of the function:

f = x1 - x2 + 3x1^2 + 2x1x2 + x2^2 + 1

However I would also like to add a constraint to this so that x1*x2 > 0

In trying to implement this I have looked at the book engineering optimization by S.S. Rao and one method shown has sort to use a penalty function 'r', so that the function now equals: f(x1,x2) - r*(1/g(x1,x2)).

To this end i have recalculated gradf to give:

gradf = [1 + 6*x(1) + x(2)*(2 - r/((x(1)*x(2)-1)^2) ); -2 + 2*x(2) + x(1)*(2 - r /((x(1)*x(2) -1 )^2)];

and tried values of r ranging from 0 ---> 1000

However I still end up with a minimum value inside the constraint.

Any kind person out there willing to lend some assistance?

% minimise f = x1 - x2 + 3x1^2 + 2x1x2 + x2^2 + 1

clear

clf

x = [2 1]'; % initial solution

xx = x';

ll = [];

ss = [];

for icount = 1:30

gradf = [1 + 6*x(1) + 2*x(2) ; -2 + 2*x(1) + 2*x(2)];

% gradf = [1 + 4*x(1) + 2*x(2); -1 + 2*x(1) + 2*x(2)]

s = -gradf;

% minimise f(x + lamda*s)

a = x(1); b = s(1); c= x(2); d = s(2);

% Alambda^2 + Blambda + C

C = a - 2*c + 3*a^2 + 2*a*c + c^2 + 1;

B = b - 2*d + 6*a*b + 2*a*d + 2*b*c + 2*c*d;

A = 3*b^2 + 2*b*d + d^2;

lambda = -B/(2*A); % optimum step length

x = x + lambda * s;

xx = [xx;x'];

ll = [ll;lambda];

ss = [ss;s'];

end

xx = xx(1:end-1,: );

[xx ll ss]

icount = 0;

for ii =-2.5:.01:2.5

icount = icount + 1;

x1(icount) = ii;

jcount = 0;

for jj = -1.5:.01:3.5

jcount = jcount + 1;

x2(jcount) = jj;

f(icount,jcount)= ii - 2*jj + 3*ii^2 + 2*ii*jj + jj^2 + 1;

end

end

figure(1)

mesh(x2,x1,f)

figure(2)

contour(x2,x1,f)

xlabel('x2')

ylabel('x1')

hold on

plot(xx(:,2),xx(:,1),'ok-')

hold off