CODEWITHSUNDEEP
pythonnumpyscipymathematical-optimization
FooBar
FooBar

Reputation: 16508

fmin_slsqp returns different solutions to the same system

The following is as much I could boil it down.

I'm trying to solve a system of equations with 18 equations and 18 variables. For the moment, I hold 4 of these variables fixed. Originally, I got weird results. So I simplified the problem so far such that the first 9 and the last 9 equations are separate and identical. Moreover, the problem is exactly identified: There should be one unique solution.

The solution vector contains 14 elements (18 minus the 4 fixed variables). Since these are ordered properly, the first 7 solution variables should be identical to the last 7 solution variables. However, they are not.

Following is the output that I get: 

Optimization terminated successfully.    (Exit mode 0)
            Current function value: 0.125393271845
            Iterations: 18
            Function evaluations: 297
            Gradient evaluations: 18
Out[223]: 
          J         W         w         v         U         u         Y
0  0.663134  0.237578  0.251245  10.00126  0.165647  0.093939  0.906657
1  0.022635  0.825547  1.000000  10.00340  0.512898  0.089790  0.909918

Where I have stacked the first 7 variables into the first row and the next 7 variables into the second row. Clearly, these are not identical.

Code for reproduction follows

import numpy as np

# parameters

class Parameters(object):
    r = 1.03
    sBar = 0.1
    sB = 0.1
    c = 0.1
    z = 0.001
    H = 1
    epsilon = 1
    beta = 0.1

def q(theta):
    if theta <= 0:
        return 999
    return float(1)/theta

def f(theta):
    if theta < 1:
        return 0
    return 1 - float(1)/theta

# sum_all allows to see residual as vector, not summed
def twoSectorFake(x, Param, sum_all=True):
    JBar, WBar, wBar, vBar, UBar, uBar, YBar = x[:7]
    JB, WB, wB, vB, UB, uB, YB = x[7:]
    VBar = 0
    VB = 0
    pB = 1
    pBar = 1
    #theta = float(vB + vBar)/u
    thetaBar = float(vBar)/uBar
    thetaB = float(vB)/uB
    res = np.empty(18,)
    res[0] = Param.r*JBar - (pBar - wBar - Param.sBar*(JBar - VBar) )
    res[1] = Param.r * VBar - ( -Param.c + q(thetaBar) * (JBar - VBar) )
    res[2] = Param.r * WBar - (wBar - Param.sBar * (WBar - UBar) )
    res[3] = Param.r * UBar - (Param.z + f(thetaBar) * (WBar - UBar) )
    res[4] = Param.sBar * YBar - vBar * q(thetaBar)
    res[5] = Param.sBar * YBar - uBar * f(thetaBar)
    res[6] = JBar - (1 - Param.beta) * (JBar + WBar - UBar)
    res[7] = Param.H - YBar - uBar
    res[8] = thetaBar * uBar - vBar

    res[9] = Param.r*JB - (pB - wB - Param.sB*(JB - VB))
    res[10] = Param.r*VB - ( -Param.c + q(thetaB) * (JB - VB))
    res[11] = Param.r * WB - (wB - Param.sB * (WB - UB))
    res[12] = Param.r * UB - (Param.z + f(thetaB) * (WB - UB))
    res[13] = Param.sB * YB - vB*q(thetaB)
    res[14] = Param.sB * YB - uB * f(thetaB)
    res[15] = JB - (1 - Param.beta) * (JB + WB - UB)
    res[16] = Param.H - YB - uB
    res[17] = thetaB * uB - vB

    idx = abs(res > 10000)
    # don't square too big numbers, they may become INF and the problem won't solve
    res[idx] = abs(res[idx])
    res[~idx] = res[~idx]**2
    if (sum_all==False):
        return res
    return sum(res)

Param = Parameters()

x2 = np.empty(0,)
boundaries2 = []
# JBar
x2 = np.append(x2, 1)
boundaries2.append([0, 100])
# WBar
x2 = np.append(x2, 1)
boundaries2.append([0, 100])
# wBar
x2 = np.append(x2, 0.5)
boundaries2.append([0.01, 100])
# vBar
x2 = np.append(x2, 10)
boundaries2.append([0.01, 100000])
# UBar
x2 = np.append(x2, float(Param.z)/(Param.r-1)+1)
boundaries2.append([float(Param.z)/(Param.r-1) - 0.1, 100])
# uBar
x2 = np.append(x2, 0.5*Param.H)
boundaries2.append([0.0000001, Param.H])
# YBar
x2 = np.append(x2, 0.5*Param.H)
boundaries2.append([0.0001, Param.H])
# JB
x2 = np.append(x2, 1)
boundaries2.append([0, 100])
# WB
x2 = np.append(x2, 1)
boundaries2.append([0, 100])
# wB
x2 = np.append(x2, 0.5)
boundaries2.append([1, 100])
# vB
x2 = np.append(x2, 10)
boundaries2.append([0.01, 100])
# UB
x2 = np.append(x2, float(Param.z)/(Param.r-1)+1)
boundaries2.append([float(Param.z)/(Param.r-1) - 0.1, 100])
# uB
x2 = np.append(x2, 0.5*Param.H)
boundaries2.append([0.0000001, Param.H])
# YB
x2 = np.append(x2, 0.5*Param.H)
boundaries2.append([0.0001, Param.H])

result = optimize.fmin_slsqp(func=twoSectorFake, x0=x2, bounds=boundaries2, args=(Param,), iter=200)
res1 = result[:7]
res2 = result[7:]
df = pd.DataFrame(np.reshape(res1, (1,7)), columns=['J', 'W', 'w', 'v', 'U', 'u','Y'])
df.loc[1] = np.reshape(res2, (1,7))
print df

Upvotes: 2

Views: 112

Answers (1)

gg349
gg349

Reputation: 22701

EDIT:

your variable boundaries2 isn't symmetric, i.e. boundaries2[7:]!=boundaries2[:7] .

Try writing 

boundaries2 = boundaries2[7:]*2

just before your call to fmin_slsqp and you get a symmetric local minimum. I leave the previous general comments on your setup below, because they apply in any case.

First of all, are you sure that there exists only one solution to your problem? if not, I wouldn't expect necessarily the scipy numerical routine to return the solution you have in mind. It can converge to any other non-symmetric solution.

In second instance, you are not solving the system of equations. If you evaluate twoSectorFake(result,Param) at the end of your code, you get 0.15. You may be better off with other root solvers, see the root finding section.

This means that you're then looking at a local minimum of your target function, i.e. not a zero. Again, there's no reason why the numerical routine should necessarily calculate a symmetric local minimum.

Upvotes: 1

Related Questions