我有一个非线性方程组。我不能很好地猜测初值。我希望至少有一组正根因为在经济学中，这些变量的负值没有多大意义。
# -*- coding: utf-8 -*-
"""
Created on Sat Oct 15 21:48:56 2016
@author: Nick
"""
import scipy as sp
from scipy.optimize import root, fsolve
import numpy as np
#from scipy.optimize import *
el = 1.1
eg = el
ej = 10
om = 0.3
omg = 0.3
rhog = 0.8
xi = 0.9
mun = 2
pidss = 0.02
muc = 0.001
ec = 2.00 # sims obtains 2.47
beta = 0.998
h = 0.8
kappa = 4.00
n = 1/3.0
alpha = 1/3.0
delta = 0.025
egs = eg
oms = 0.2
omgs = oms
rhom = 0.7
psiygap = 1.000
psipi = 2.500
rhoicu = 0.800
taudss = 0.01 # steady state tax on domestic consumption (setting it as 0 would create algebraic difficulties)
taumss = 0.01 # steady state tax on imported consumption for domestic country
taukss = 0.01 # steady state tax on rental income from capital for domestic country block
taunss = 0.01 # steady state tax on labor for domestic country
tauydss = 0.05
gss = 0.23 # steady state government spending as a propostion of gdp for domestic country block
gsss = 0.23 # steady state government spending as a propostion of gdp for foreign country block
taudsss = 0.01
taumsss = 0.01
tauksss = 0.01
taunsss = 0.01
tauydsss = 0.01 # steady state tax rate on output for foreign country block
tauss = 1.0 # Steady state terms of trade
icu = ((1+pidss)/beta) - 1
mc = ((ej - 1)/ej)
r = (1/taukss) * ((1/beta) - (1-delta))
rs = (1-tauksss) * ((1/beta) - (1-delta))
KN = (mc*alpha/r)**(1/(1-alpha))
KNs = (mc*alpha/rs)**(1/(1-alpha))
psigma = (1-xi) * (1/(1-tauydss) - xi)**(-1)
psigmas = (1-xi) * (1/(1-tauydsss) - xi)**(-1)
w = (1-alpha) * mc * (KN)**(alpha)
z = np.zeros(16)
def fun(z):
Yd = z[0]
N = z[1]
X = z[2]
I = z[3]
Cd = z[4]
Cm = z[5]
Gd = z[6]
Gm = z[7]
Yds = z[8]
Ns = z[9]
Xs = z[10]
Is = z[11]
Cds = z[12]
Cms = z[13]
Gds = z[14]
Gms = z[15]
print (z)
f = np.zeros(16)
f[0] = N - ( (X - muc)**(-ec) * ((1-alpha)/(mun)) * (mc)**(1/(1-alpha)) * (alpha/r)** (1-taunss) )
f[1] = Yd - ( Cd + Gd + I + ((1-n)/n) *(Cms + Gms) )
f[2] = Yd - ( (KN)**(alpha) * (psigma/(1-tauydss)**(ej)) )
f[3] = Cd - ( X * ((1-om)/(1+taudss)**(el)) *((1-om)*(1+taudss)**(1-el) + om * (1+taumss)**(1-el) * tauss**(1-el))**(el/(1-el)) )
f[4] = Gd - ( ((gss*Yd * (1-omg))/(1+taudss)**(eg) ) *((1-omg)*(1+taudss)**(1-eg) + omg* (1+taumss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg) ) )
f[5] = I - ( delta* KN * N )
f[6] = Cm -( (X * (1-om)/(1+tauydss)**(el) ) *((1-om)*(1+taudss)**(1-el) + om* (1+taumss)**(1-el) * tauss**(1-el))**(el/(1-el)) )
f[7] = Gm - ( ((gss*Yd * (omg))/(1+taumss)**(eg) ) *((1-omg)*(1+taudss)**(1-eg) + omg* (1+taumss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg) ) )
f[8] = Ns - ( (Xs - muc)**(-ec) * ((1-alpha)/(mun)) * (mc)**(1/(1-alpha)) * (alpha/rs)** (1-taunsss) )
f[9] = Yds - ( Cds + Gds + Is + (n/(1-n)) *(Cm + Gm) )
f[10] = Yds - ( (KNs)**(alpha) * (psigmas/(1-tauydsss)**(ej) ) )
f[11] = Cds - ( Xs * ((1-oms)/(1+taudsss)**(el))* ((1-oms)*(1+taudsss)**(1-el) + oms* (1+taumsss)**(1-el) * tauss**(1-el))**(el/(1-el)) )
f[12] = Gds - ( ((gsss*Yds * (1-omgs))/(1+taudsss)**(eg) ) *((1-omgs)*(1+taudsss)**(1-eg) + omgs* (1+taumsss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg) ) )
f[13] = Is - ( delta* KNs * Ns )
f[14] = Cms -( (Xs * (1-oms)/(1+tauydsss)**(el) ) *((1-oms)*(1+taudsss)**(1-el) + oms* (1+taumsss)**(1-el) * tauss**(1-el))**(el/(1-el)) )
f[15] = Gms - ( ((gsss*Yds * (omgs))/(1+taumsss)**(eg) ) *((1-omgs)*(1+taudsss)**(1-eg) + omgs* (1+taumsss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg) ) )
return f
z = sp.optimize.root(fun, [100,100,70,30,50,20,50,20,100,100,100,100,100,100,100,100], method='lm')
#z = fsolve(fun, [0,0,0.0,0,1,1,1,1,1,1,1,1,1,1,1,1])
print(z)
解为以下根
success: True
x: array([ 3.64725445e-01, 1.02848541e-06, -1.86761721e+02,
9.52089296e-10, -1.30733205e+02, -1.25265418e+02,
5.87207967e-02, 2.51660557e-02, 3.36422990e+00,
5.18324506e-04, 8.17060628e+01, 4.87111630e-04,
6.53648502e+01, 6.53648502e+01, 6.19018302e-01,
1.54754576e-01])
给定根的初始估计值，数值寻根算法在变量空间中沿一定方向移动，直到找到根为止。显然，使用这种方法，没有必要要求返回的根在一定的区间内有界——这完全取决于初始估计的好坏(以及算法使用的搜索方法)。
另一种可以返回有界根的方法是将寻根问题作为优化(例如最小化)问题，因为在优化问题中提供约束是有意义的。但是，您必须提供一个适当的目标函数，其最小值出现在原始函数的根处(对于这样的函数有许多选项，通常选择是启发式的)。
一个这样的函数是平方和f[0]**2 + f[1]**2 + ... + f[15]**2。显然，这个函数的最小值为0，当和的每一个单独的项都为0时，即在它们的根处。您可以使用Scipy的least_squares来执行这种最小化，它还允许为优化变量提供边界。
在变量没有任何边界的情况下，使用相同的初始根估计，least_squares返回与root相同的解:
from scipy.optimize import least_squares
z_ls = least_squares(fun, [100,100,70,30,50,20,50,20,100,100,100,100,100,100,100,100])
print(z_ls.x)
print(z_ls.cost)
[ 3.6473e-01 1.0285e-06 -1.8676e+02 9.5209e-10 -1.3073e+02 -1.2527e+02 5.8721e-02 2.5166e-02 3.3642e+00 5.1832e-04 8.1706e+01 4.8711e-04 6.5365e+01 6.5365e+01 6.1902e-01 1.5475e-01] 4.16527754459e-26
(注意z_ls.cost是此时的平方和，在数值精度范围内为0)
现在，使用least_squares约束估计是非负的:
z_ls = least_squares(fun, [100,100,70,30,50,20,50,20,100,100,100,100,100,100,100,100],bounds = (0,np.inf))
print(z_ls.x)
print(z_ls.cost)
[ 5.9581e-01 2.1229e+01 4.2108e-02 1.1820e-37 2.0493e-33 1.1914e-33 5.8857e-37 9.3812e-37 3.4508e+00 2.5054e+00 1.1516e+00 2.2395e+00 8.0630e-01 4.2258e-01 5.1994e-01 1.3867e-37] 0.237262813475
返回的估计值确实包含非负元素。然而，z_ls.cost(显著)大于零，表明该解是而不是根。这意味着两者之一:
- 初始点不足以导致非负根。
- 这个问题不存在非负根。
如果您对上述没有任何见解，您唯一可以做的就是尝试不同的初始化值，并希望返回所需的根(直接通过root或像上面那样的最小化公式)。
假设您想找到始终为正的变量v的值。
让求解器找到它的对数(log_v = log(v))，然后设置v = exp(log_v)，这样当求解器在整个域中搜索log_v时，v总是正的。
def fun(log_v):
v = exp(log_v)
...[do the same math using v]
sol = fsolve(fun, (1,))
v = exp(sol[0])