Trying to find a stationary point using the active set method

46 Views Asked by At

I would like if someone could verify if i did everything ok or not.

So i want to minimize function f using the active set method.

$$ min f(x_{1},x_{2}) = 2x_{1}^2 + 3x_{1}x_{2}^2 + 6x_{1}^2x_{2} +x_{2}^2 + 4x_{1}x_{2} + x_{2} $$

subject to

$$ x_{1} + x_{2} <= 32 \iff -x_{1} - x_{2} >= -32, constraint 1$$

$$ x_{1} + 2x_{2} >= 8 , constraint 2$$

$$ x_{1} >= 0 , constraint 3$$

$$ x_{2} >= 0 , constraint 4$$

The gradient and hessian of f are:

$$gF =\begin{bmatrix}4x_{1}+3x_{2}^2+12x_{1}x_{2}+4x_{2} \\ 6x_{1}x_{2}+6x_{1}^2+2x_{2}+4x_{1}+1 \end{bmatrix}$$

$$hF =\begin{bmatrix}4 + 12x_{2} & 6x_{2} + 12x_{1} + 4 \\ 6x_{2} + 12x_{1}+4 & 6x_{1} +2 \end{bmatrix}$$

So to begin with, in the initial point there are no active constraints, which means that the problem is locally unconstrained.

$$ p = -hF(x)^{-1} \times gF(x) $$

Now to see if that step decreases the function i calculated the hessian eigen values:

$$eig(hF(x)) =\begin{bmatrix}-28.85 \\ 88.85 \end{bmatrix}$$

It's not positive definite so i will rather use the step of steepest descent

$$p = -gF(x) =\begin{bmatrix}-159 \\ -127 \end{bmatrix}$$

Now to compute the step size and make sure the step is in the feasible zone:

With alfa1 (step size for constraint 1) i have noticed that this gives negative:

$$ -a_{1}\times p = -\begin{bmatrix}1 & 2 \end{bmatrix}\times p = -286$$

So i dont have to worry about this inactive constraint for now

For constraint 2:

$$ alfa2 = \frac{\begin{bmatrix}1 & 2 \end{bmatrix} \times x - 8}{-\begin{bmatrix}1 & 0 \end{bmatrix}\times p} = 0.0024$$

For constraint 3 and 4 i have:

$$ alfa3 = \frac{\begin{bmatrix}1 & 0 \end{bmatrix} \times x - 0}{-\begin{bmatrix}1 & 0 \end{bmatrix}\times p} = 0.0189$$

$$ alfa4 = \frac{\begin{bmatrix}0 & 1 \end{bmatrix} \times x - 0}{-\begin{bmatrix}0 & 1 \end{bmatrix}\times p} = 0.0236$$

alfa2 has the minimum step to the boundary, so:

$$ x = x + alfa2 \times p = \begin{bmatrix}2.6150 \\ 2.6925 \end{bmatrix}$$

So now constraint 2 is active.

Let A be the matrix for active constraints and Z its null matrix

$$ A = \begin{bmatrix}1 & 2 \end{bmatrix} $$ $$ Z = \begin{bmatrix}-0.8944 \\ 0.4472 \end{bmatrix} $$

I want to use the newton step but first i need to make sure if it leads to a decrease of function f.

$$eig(hF(x)) =\begin{bmatrix}-25.3693 \\ 79.3693 \end{bmatrix}$$

It's not positive definite so i will need to use the direction of steepest descent again.

Reduced Steepest-Descent Direction:

$$p = Z \times Z^{T} \times gF(x) = \begin{bmatrix}-61.9275 \\ 30.9637 \end{bmatrix}$$

Now to compute the step size and make sure the step is in the feasible zone:

For constraint 1:

$$ -\begin{bmatrix}1 & 1 \end{bmatrix}\times p = -30.9637$$

So i dont have to worry about this inactive constraint for now.

For constraint 3:

$$ alfa3 = \frac{\begin{bmatrix}1 & 0 \end{bmatrix} \times x - 0}{-\begin{bmatrix}1 & 0 \end{bmatrix}\times p} = 0.0422$$

For constraint 4:

$$ -\begin{bmatrix}0 & 1 \end{bmatrix}\times p = -30.9637$$

Dont have to worry about this inactive constraint for now.

So the step size will be 0.0422

$$ x = x + alfa3 \times p = \begin{bmatrix}0 \\ 4 \end{bmatrix}$$

Now the matrix for active constraints and its null matrix are:

$$ A = \begin{bmatrix}1 & 2 \\ 1 & 0 \end{bmatrix} $$ Z is an empty matrix

If Z is an empty matrix then the reduced gradient also is.

With that i computed the lagrange multipliers

$$ \lambda = A^{T} \times gF(x)^{-1} = \begin{bmatrix} 4.5 \\ 59.5 \end{bmatrix}$$

With this i can say i found its minimum.

Also

$$f(\begin{bmatrix} 3 \\ 3 \end{bmatrix}) = 309$$

and

$$f(\begin{bmatrix} 0 \\ 4 \end{bmatrix}) = 20 $$

I just find this problem weird in the sense that i only needed 2 iterations to solve it. I made all these computations in matlab.

If you guys wanna see it:

clc
clear

f = @(x) 2 * x(1)^2 + 3*x(1)*(x(2)^2) + 6*(x(1)^2)*x(2) + x(2)^2 + 4*x(1)*x(2) + x(2);
f_gradiente = @(x) [4*x(1)+3*(x(2)^2)+12*x(1)*x(2)+4*x(2);6*x(1)*x(2)+6*(x(1)^2)+2*x(2)+4*x(1)+1];
f_hessiana = @(x) [4+12*x(2),6*x(2)+12*x(1)+4;6*x(2)+12*x(1)+4,6*x(1)+2];
constraint_1 = [-1,-1];
constraint_2 = [1,2];
constraint_3 = [1,0];
constraint_4 = [0,1];
b1 = -32;
b2 = 8;
b3 = 0;
b4 = 0;

% Iteracao 1

% Nao ha constraints ativas
x = [3;3];

f_gradiente(x)
passo_newton = - f_hessiana(x)^(-1) * f_gradiente(x) % como nao ha constraints, esta iteracao os passos a usar nao sao reduzidos
f_hessiana(x)
eig(f_hessiana(x)) % da negativo, nao é positive definite, n se pode usar o passo de newton

% Usar antes o passo de cauchy
p = -f_gradiente(x)

denominador_alfa = - constraint_1 * p % da positivo o simetrico disto, n e necessario calcular alfa
numerador_alfa = constraint_1 * x -b1;

alfa1 = numerador_alfa / denominador_alfa

denominador_alfa = - constraint_2 * p 
numerador_alfa = constraint_2 * x -b2;

alfa2 = numerador_alfa / denominador_alfa

denominador_alfa = - constraint_3 * p
numerador_alfa = constraint_3 * x -b3;

alfa3 = numerador_alfa / denominador_alfa

denominador_alfa = - constraint_4 * p
numerador_alfa = constraint_4 * x -b4;

alfa4 = numerador_alfa / denominador_alfa

% mais pequeno e o alfa2 com 0.0024

x = x + alfa2*p

% Iteracao 2

% A constraint 2 esta ativa

A = [1,2];
Z = null(A)

f_hessiana(x)
eig(f_hessiana(x)) % da negativo

p = Z * -Z' * f_gradiente(x) % passo de cauchy reduzido

denominador_alfa = - constraint_1 * p 
numerador_alfa = constraint_1 * x -b1;

alfa1 = numerador_alfa / denominador_alfa %da positivo o simetrico disto, n e necessario calcular alfa

denominador_alfa = - constraint_2 * p 
numerador_alfa = constraint_2 * x -b2;

alfa2 = numerador_alfa / denominador_alfa %da positivo o simetrico disto, n e necessario calcular alfa, nota: numero negativo mas perto de 0

denominador_alfa = - constraint_3 * p
numerador_alfa = constraint_3 * x -b3;

alfa3 = numerador_alfa / denominador_alfa

denominador_alfa = - constraint_4 * p %da positivo o simetrico disto, n e necessario calcular alfa

% mais pequeno e o alfa3 com 0.0422

x = x + alfa3*p

% Iteracao 3

% A constraint 2 e 3 estao ativas

A = [1,2;1,0];
Z = null(A)

multiplicador_lagrange = A' \ f_gradiente(x) % tudo positivo, acaba aqui

Edit: Made a mistake in the initial active constraint matrix, had to do everything again.