CODEWITHSUNDEEP
pythonmath
Maksim Vorontsov
Maksim Vorontsov

Reputation: 937

Penalty function method

I am trying to implement penalty function method for minimizing function. I need to find the minimum of Rosenbrok's function.

I am using this penalty function:

Penalty Function

Penalty Function

First of all, I have found the minimum using scipy.optimize.minimize:

from scipy.optimize import minimize, rosen
rz = lambda x: (1-x[0])**2 + 100*(x[1] - x[0]**2)**2;
h_1 = lambda x: (x[0] - 2 * x[1] + 2);
h_2 = lambda x: (-x[0] - 2 * x[1] + 6);
h_3 = lambda x: (-x[0] + 2 * x[1] + 2);

x0 = [2.3, 5];
cons = ({'type': 'ineq', 'fun': h_1},
       {'type': 'ineq', 'fun': h_2},
       {'type': 'ineq', 'fun': h_3}) 
minimize(rz, x0, constraints=cons)

The answer is x: array([ 0.99971613, 0.99942073])

Then I am trying to find the minimum using my implementation of penalty method:

x_c = [2.3, 3];
i = 1;
while i < 1000:
    curr_func = lambda x: rz(x) + i*(h_1(x)**2 + h_2(x)**2 + h_3(x)**2)
    x_c = minimize(curr_func, x_c).x;
    i  *= 1.2;
print(answer.x);

Which gives me [ 2.27402022 1.4157964 ] (if I increase the number of iterations, final values are even greater).

Where is the mistake in my implementation? Thanks.

P.S. Function curr_func is specific for my constraints, of course, when they are all 'inequals' type.

Upvotes: 4

Views: 6777

Answers (1)

arghbleargh
arghbleargh

Reputation: 3170

The problem you have is that the h_i in your formula are for equality constraints, whereas the problem you are solving is for inequality constraints, which correspond to the g_i in your formula. Hence, your penalty function should be using terms like min(0, h_1(x))**2 instead of h_1(x)**2. To see why this is the case, just think about what happens if i = 1000 and x is the desired solution (1, 1). Then, the penalty will include a term i * h_1(x)**2 = 1000, which is huge.

Note that I used min instead of max because it seems like the inequality you want to enforce is h_1(x) >= 0. That means as long as h_1(x) >= 0, the penalty should be zero, but as soon as h_1(x) goes negative, you start penalizing. If it's actually h_1(x) <= 0 you want, then you use max (then you'll have to switch h_1 with -h_1 when you use scipy.optimize.minimize).

BTW, since i is usually an index variable, it's probably better to name the penalty weight something else, like a.

Upvotes: 5

Related Questions