将 GEKKO 与快速傅里叶变换结合使用
Using GEKKO with Fast Fourier Transform
我实际上正在尝试使用 GEKKO 上可用的 IPOPT 优化器来优化大型非凸和非线性 problem.In 为了做到这一点我需要使用快速傅立叶变换 scipy.First，让我们修复示例数据（为简单起见）：
import numpy as np
import pandas as pd
import math
from scipy import *
from scipy.integrate import quad
import scipy.stats as ss
import scipy.optimize as scpo
from scipy import sparse
from scipy.fftpack import ifft
from scipy.interpolate import interp1d
from scipy.optimize import fsolve
from functools import partial
r,prices, strikes, spreads,s0,T = 0,array([1532.45 , 1507.55 , 1482.65 , 1457.8 , 1432.95 , 1408.1 ,
1383.25 , 1358.45 , 1333.6 , 1308.8 , 1284. , 1259.2 ,
1234.45 , 1209.7 , 1037.15 , 1012.55 , 988.05 , 963.35 ,
938.9 , 914.3 , 889.8 , 865.5 , 841. , 816.65 ,
792.45 , 768.1 , 743.95 , 719.85 , 695.85 , 672. ,
648.1 , 624.5 , 600.9 , 577.5 , 554.2 , 531. ,
508.15 , 485.35 , 462.9 , 440.65 , 418.65 , 396.85 ,
375.55 , 354.5 , 333.9 , 313.65 , 293.85 , 255.75 ,
237.55 , 219.8 , 202.8 , 186.35 , 170.55 , 155.4 ,
141.05 , 127.4 , 113.5 , 101.35 , 90.1 , 79.65 ,
70. , 61.3 , 53.4 , 46.35 , 34.5 , 29.6 ,
25.35 , 18.55 , 15.85 , 13.55 , 11.55 , 9.9 ,
7.35 , 5.45 , 3.075, 2.7 ]),array([12500., 12525., 12550., 12575., 12600., 12625., 12650., 12675.,
12700., 12725., 12750., 12775., 12800., 12825., 13000., 13025.,
13050., 13075., 13100., 13125., 13150., 13175., 13200., 13225.,
13250., 13275., 13300., 13325., 13350., 13375., 13400., 13425.,
13450., 13475., 13500., 13525., 13550., 13575., 13600., 13625.,
13650., 13675., 13700., 13725., 13750., 13775., 13800., 13850.,
13875., 13900., 13925., 13950., 13975., 14000., 14025., 14050.,
14075., 14100., 14125., 14150., 14175., 14200., 14225., 14250.,
14300., 14325., 14350., 14400., 14425., 14450., 14475., 14500.,
14550., 14600., 14700., 14725.]),array([29.7 , 29.7 , 29.7 , 29.8 , 29.7 , 29.8 , 29.7 , 29.7 , 29.8 ,
29.8 , 29.8 , 29.8 , 29.7 , 29.8 , 10.3 , 10.3 , 10.5 , 10.3 ,
10.6 , 10.4 , 10.4 , 10.6 , 10.4 , 10.5 , 10.7 , 10.4 , 10.5 ,
10.5 , 10.5 , 10.8 , 10.6 , 10.8 , 10.6 , 10.8 , 10.6 , 10.6 ,
10.9 , 10.5 , 10.8 , 10.7 , 10.7 , 10.3 , 10.5 , 10.4 , 10.2 ,
10.1 , 9.9 , 9.5 , 9.3 , 9. , 8.8 , 8.5 , 8.3 , 8.2 ,
7.9 , 7.6 , 3.8 , 3.7 , 3.6 , 3.5 , 3.4 , 3.2 , 3.2 ,
3.1 , 2.8 , 2.6 , 2.5 , 2.3 , 2.1 , 2.1 , 1.9 , 1.8 ,
1.7 , 1.5 , 1.25, 1.2 ]),14000,0.05
然后傅里叶函数：
class Heston_pricer():
def __init__(self, Option_info, Process_info ):
"""
Process_info: a instance of "Heston_process.", which contains (mu, rho, sigma, theta, kappa)
Option_info: of type Option_param, which contains (S0,K,T)
"""
self.r = Process_info.mu # interest rate
self.sigma = Process_info.sigma # Heston parameters
self.theta = Process_info.theta
self.kappa = Process_info.kappa
self.rho = Process_info.rho
self.S0 = Option_info.S0 # current price
self.v0 = Option_info.v0 # spot variance
self.K = Option_info.K # strike
self.T = Option_info.T # maturity(in years)
self.exercise = Option_info.exercise
self.payoff = Option_info.payoff
# payoff function
def payoff_f(self, S):
if self.payoff == "call":
Payoff = np.maximum( S - self.K, 0 )
elif self.payoff == "put":
Payoff = np.maximum( self.K - S, 0 )
return Payoff
# FFT method. It returns a vector of prices.
def FFT(self, K): # K: strikes
K = np.array(K)
# Heston characteristic function (proposed by Schoutens 2004)
def cf_Heston_good(u, t, v0, mu, kappa, theta, sigma, rho):
xi = kappa - sigma*rho*u*1j
d = m.sqrt( xi**2 + sigma**2 * (u**2 + 1j*u) )
g1 = (xi+d)/(xi-d)
g2 = 1/g1
cf = m.exp( 1j*u*mu*t + (kappa*theta)/(sigma**2) * ( (xi-d)*t - 2*m.log( (1-g2*m.exp(-d*t))/(1-g2) ))\
+ (v0/sigma**2)*(xi-d) * (1-m.exp(-d*t))/(1-g2*m.exp(-d*t)))
return cf
cf_H_b_good = partial(cf_Heston_good, t=self.T, v0=self.v0, mu=self.r, theta=self.theta,
sigma=self.sigma, kappa=self.kappa, rho=self.rho)
if self.payoff == "call":
return fft_(K, self.S0, self.r, self.T, cf_H_b_good)
elif self.payoff == "put": # put-call parity
return fft_(K, self.S0, self.r, self.T, cf_H_b_good) - self.S0 + K*m.exp(-self.r*self.T)
class Heston_process():
def __init__(self, mu=0, rho=0, sigma=0.00001, theta=0.4, kappa=.00001):
"""
r: risk free constant rate
rho: correlation between stock noise and variance noise (|rho| must be <=1)
theta: long term mean of the variance process(positive)
sigma: volatility coefficient(positive)
kappa: mean reversion coefficient for the variance process(positive)
"""
self.mu, self.rho, self.theta, self.sigma, self.kappa = mu, rho, theta, sigma, kappa
def fft_(K, S0, r, T, cf): # interp support cubic
"""
K = vector of strike
S0 = spot price scalar
cf = characteristic function
"""
N=2**15 # FFT more efficient for N power of 2
B = 500 # integration limit
dx = B/N
x = np.arange(N) * dx
weight = 3 + (-1)**(np.arange(N)+1) # Simpson weights
weight[0] = 1; weight[N-1]=1
dk = 2*np.pi/B
b = N * dk /2
ks = -b + dk * np.arange(N)
integrand = m.exp(- 1j * b *(np.arange(N))*dx) * cf(x - 0.5j) * 1/(x**2 + 0.25) * weight * dx/3
integral_value = np.real(ifft(integrand)*N)
spline_cub = interp1d(ks, integral_value, kind="cubic") # cubic will fit better than linear
prices = S0 - m.sqrt(S0 * K) * m.exp(-r*T)/np.pi * spline_cub( m.log(S0/K) )
return prices
# A class that stores option parameters (in order to write BS/Heston class neatly)
class Option_param():
def __init__(self, S0=10000, K=10000, T=.1, v0=0.04, payoff="call", exercise="European"):
"""
S0: current stock price
K: Strike price
T: time to maturity
v0: (optional) spot variance
exercise: European or American
"""
self.S0, self.v0, self.K, self.T, self.exercise, self.payoff = S0, v0, K, T, exercise, payoff
现在让我们使用 GEKKO：
#Initialize Model
m = GEKKO()
#define parameter
eq = m.Param(value=5)
#initialize variables
x1,x2,x3,x4,x5 = [m.Var(lb=-1, ub=1),m.Var(lb=1e-3, ub=1),m.Var(lb=1e-3, ub=1),m.Var(lb=1e-3, ub=20),m.Var(lb=1e-3, ub=1)]
#initial values
x1.value = 0
x2.value = 0.5
x3.value = 0.5
x4.value = 0.5
x5.value = 0.5
X = [x1,x2,x3,x4,x5]
#Equations,Feller Condition
m.Equation(2*x3*x4 - x2*x2 >=0)
def least_sq(x):
""" Objective function """
Heston_param = Heston_process(mu=0, rho=X[0], sigma=X[1], theta=X[2], kappa=X[3])
m = 1/(spreads**2)
if len(m) == 1:
l = 1
else:
l = (m - np.min(m))/(np.max(m)-np.min(m))
opt_param = Option_param(S0=s0, K=strikes, T=T, v0=X[4], exercise="European", payoff="call" )
Hest = Heston_pricer(opt_param, Heston_param)
prices_calib = Hest.FFT(strikes)
results = (l * (prices_calib-prices)**2)/len(prices)
return results
m.Obj(m.sum(least_sq(X)))
#Set global options
m.options.IMODE = 3 #steady state optimization
#Solve simulation
m.solve()
#Results
print('')
print('Results')
print('x1: ' + str(x1.value))
print('x2: ' + str(x2.value))
print('x3: ' + str(x3.value))
print('x4: ' + str(x4.value))
print('x5: ' + str(x5.value))
这里的问题是 ifft scipy 函数不起作用，因为 GEKKO.The 给出的变量类型不同问题是我不知道如何替换或避免它.
错误是这样的：
<ipython-input-16-305a2bb0769b> in fft_(K, S0, r, T, cf)
78
79 integrand = m.exp(- 1j * b *m.CV(np.arange(N))*dx) * cf(x - 0.5j) * 1/(x**2 + 0.25) * weight * dx/3
---> 80 integral_value = np.real(ifft(integrand)*N)
81 spline_cub = interp1d(ks, integral_value, kind="cubic") # cubic will fit better than linear
82 prices = S0 - m.sqrt(S0 * K) * m.exp(-r*T)/np.pi * spline_cub( m.log(S0/K) )
/usr/local/lib/python3.7/dist-packages/scipy/_lib/deprecation.py in call(*args, **kwargs)
18 warnings.warn(msg, category=DeprecationWarning,
19 stacklevel=stacklevel)
---> 20 return fun(*args, **kwargs)
21 call.__doc__ = msg
22 return call
<__array_function__ internals> in ifft(*args, **kwargs)
/usr/local/lib/python3.7/dist-packages/numpy/fft/_pocketfft.py in ifft(a, n, axis, norm)
274 a = asarray(a)
275 if n is None:
--> 276 n = a.shape[axis]
277 if norm is not None and _unitary(norm):
278 inv_norm = sqrt(max(n, 1))
IndexError: tuple index out of range
谁能帮我调试this.Thank你
Gekko 要求表达式不是黑盒，而是可以用特殊类型的变量（Gekko 类型）来表示，用于自动微分和稀疏检测。使用 Scipy.optimize.minimize 等求解器可能会更好地解决这个问题。这是一个comparison of the two on a simple problem.
Scipy
import numpy as np
from scipy.optimize import minimize
def objective(x):
return x[0]*x[3]*(x[0]+x[1]+x[2])+x[2]
def constraint1(x):
return x[0]*x[1]*x[2]*x[3]-25.0
def constraint2(x):
sum_eq = 40.0
for i in range(4):
sum_eq = sum_eq - x[i]**2
return sum_eq
# initial guesses
n = 4
x0 = np.zeros(n)
x0[0] = 1.0
x0[1] = 5.0
x0[2] = 5.0
x0[3] = 1.0
# show initial objective
print('Initial Objective: ' + str(objective(x0)))
# optimize
b = (1.0,5.0)
bnds = (b, b, b, b)
con1 = {'type': 'ineq', 'fun': constraint1}
con2 = {'type': 'eq', 'fun': constraint2}
cons = ([con1,con2])
solution = minimize(objective,x0,method='SLSQP',\
bounds=bnds,constraints=cons)
x = solution.x
# show final objective
print('Final Objective: ' + str(objective(x)))
# print solution
print('Solution')
print('x1 = ' + str(x[0]))
print('x2 = ' + str(x[1]))
print('x3 = ' + str(x[2]))
print('x4 = ' + str(x[3]))
盖柯
from gekko import GEKKO
m = GEKKO()
#initialize variables
x1,x2,x3,x4 = [m.Var(lb=1,ub=5) for i in range(4)]
x1.value = 1; x2.value = 5; x3.value = 5; x4.value = 1
m.Equation(x1*x2*x3*x4>=25)
m.Equation(x1**2+x2**2+x3**2+x4**2==40)
m.Minimize(x1*x4*(x1+x2+x3)+x3)
m.solve()
print('x1: ' + str(x1.value))
print('x2: ' + str(x2.value))
print('x3: ' + str(x3.value))
print('x4: ' + str(x4.value))
我实际上正在尝试使用 GEKKO 上可用的 IPOPT 优化器来优化大型非凸和非线性 problem.In 为了做到这一点我需要使用快速傅立叶变换 scipy.First，让我们修复示例数据（为简单起见）：
import numpy as np
import pandas as pd
import math
from scipy import *
from scipy.integrate import quad
import scipy.stats as ss
import scipy.optimize as scpo
from scipy import sparse
from scipy.fftpack import ifft
from scipy.interpolate import interp1d
from scipy.optimize import fsolve
from functools import partial
r,prices, strikes, spreads,s0,T = 0,array([1532.45 , 1507.55 , 1482.65 , 1457.8 , 1432.95 , 1408.1 ,
1383.25 , 1358.45 , 1333.6 , 1308.8 , 1284. , 1259.2 ,
1234.45 , 1209.7 , 1037.15 , 1012.55 , 988.05 , 963.35 ,
938.9 , 914.3 , 889.8 , 865.5 , 841. , 816.65 ,
792.45 , 768.1 , 743.95 , 719.85 , 695.85 , 672. ,
648.1 , 624.5 , 600.9 , 577.5 , 554.2 , 531. ,
508.15 , 485.35 , 462.9 , 440.65 , 418.65 , 396.85 ,
375.55 , 354.5 , 333.9 , 313.65 , 293.85 , 255.75 ,
237.55 , 219.8 , 202.8 , 186.35 , 170.55 , 155.4 ,
141.05 , 127.4 , 113.5 , 101.35 , 90.1 , 79.65 ,
70. , 61.3 , 53.4 , 46.35 , 34.5 , 29.6 ,
25.35 , 18.55 , 15.85 , 13.55 , 11.55 , 9.9 ,
7.35 , 5.45 , 3.075, 2.7 ]),array([12500., 12525., 12550., 12575., 12600., 12625., 12650., 12675.,
12700., 12725., 12750., 12775., 12800., 12825., 13000., 13025.,
13050., 13075., 13100., 13125., 13150., 13175., 13200., 13225.,
13250., 13275., 13300., 13325., 13350., 13375., 13400., 13425.,
13450., 13475., 13500., 13525., 13550., 13575., 13600., 13625.,
13650., 13675., 13700., 13725., 13750., 13775., 13800., 13850.,
13875., 13900., 13925., 13950., 13975., 14000., 14025., 14050.,
14075., 14100., 14125., 14150., 14175., 14200., 14225., 14250.,
14300., 14325., 14350., 14400., 14425., 14450., 14475., 14500.,
14550., 14600., 14700., 14725.]),array([29.7 , 29.7 , 29.7 , 29.8 , 29.7 , 29.8 , 29.7 , 29.7 , 29.8 ,
29.8 , 29.8 , 29.8 , 29.7 , 29.8 , 10.3 , 10.3 , 10.5 , 10.3 ,
10.6 , 10.4 , 10.4 , 10.6 , 10.4 , 10.5 , 10.7 , 10.4 , 10.5 ,
10.5 , 10.5 , 10.8 , 10.6 , 10.8 , 10.6 , 10.8 , 10.6 , 10.6 ,
10.9 , 10.5 , 10.8 , 10.7 , 10.7 , 10.3 , 10.5 , 10.4 , 10.2 ,
10.1 , 9.9 , 9.5 , 9.3 , 9. , 8.8 , 8.5 , 8.3 , 8.2 ,
7.9 , 7.6 , 3.8 , 3.7 , 3.6 , 3.5 , 3.4 , 3.2 , 3.2 ,
3.1 , 2.8 , 2.6 , 2.5 , 2.3 , 2.1 , 2.1 , 1.9 , 1.8 ,
1.7 , 1.5 , 1.25, 1.2 ]),14000,0.05
然后傅里叶函数：
class Heston_pricer():
def __init__(self, Option_info, Process_info ):
"""
Process_info: a instance of "Heston_process.", which contains (mu, rho, sigma, theta, kappa)
Option_info: of type Option_param, which contains (S0,K,T)
"""
self.r = Process_info.mu # interest rate
self.sigma = Process_info.sigma # Heston parameters
self.theta = Process_info.theta
self.kappa = Process_info.kappa
self.rho = Process_info.rho
self.S0 = Option_info.S0 # current price
self.v0 = Option_info.v0 # spot variance
self.K = Option_info.K # strike
self.T = Option_info.T # maturity(in years)
self.exercise = Option_info.exercise
self.payoff = Option_info.payoff
# payoff function
def payoff_f(self, S):
if self.payoff == "call":
Payoff = np.maximum( S - self.K, 0 )
elif self.payoff == "put":
Payoff = np.maximum( self.K - S, 0 )
return Payoff
# FFT method. It returns a vector of prices.
def FFT(self, K): # K: strikes
K = np.array(K)
# Heston characteristic function (proposed by Schoutens 2004)
def cf_Heston_good(u, t, v0, mu, kappa, theta, sigma, rho):
xi = kappa - sigma*rho*u*1j
d = m.sqrt( xi**2 + sigma**2 * (u**2 + 1j*u) )
g1 = (xi+d)/(xi-d)
g2 = 1/g1
cf = m.exp( 1j*u*mu*t + (kappa*theta)/(sigma**2) * ( (xi-d)*t - 2*m.log( (1-g2*m.exp(-d*t))/(1-g2) ))\
+ (v0/sigma**2)*(xi-d) * (1-m.exp(-d*t))/(1-g2*m.exp(-d*t)))
return cf
cf_H_b_good = partial(cf_Heston_good, t=self.T, v0=self.v0, mu=self.r, theta=self.theta,
sigma=self.sigma, kappa=self.kappa, rho=self.rho)
if self.payoff == "call":
return fft_(K, self.S0, self.r, self.T, cf_H_b_good)
elif self.payoff == "put": # put-call parity
return fft_(K, self.S0, self.r, self.T, cf_H_b_good) - self.S0 + K*m.exp(-self.r*self.T)
class Heston_process():
def __init__(self, mu=0, rho=0, sigma=0.00001, theta=0.4, kappa=.00001):
"""
r: risk free constant rate
rho: correlation between stock noise and variance noise (|rho| must be <=1)
theta: long term mean of the variance process(positive)
sigma: volatility coefficient(positive)
kappa: mean reversion coefficient for the variance process(positive)
"""
self.mu, self.rho, self.theta, self.sigma, self.kappa = mu, rho, theta, sigma, kappa
def fft_(K, S0, r, T, cf): # interp support cubic
"""
K = vector of strike
S0 = spot price scalar
cf = characteristic function
"""
N=2**15 # FFT more efficient for N power of 2
B = 500 # integration limit
dx = B/N
x = np.arange(N) * dx
weight = 3 + (-1)**(np.arange(N)+1) # Simpson weights
weight[0] = 1; weight[N-1]=1
dk = 2*np.pi/B
b = N * dk /2
ks = -b + dk * np.arange(N)
integrand = m.exp(- 1j * b *(np.arange(N))*dx) * cf(x - 0.5j) * 1/(x**2 + 0.25) * weight * dx/3
integral_value = np.real(ifft(integrand)*N)
spline_cub = interp1d(ks, integral_value, kind="cubic") # cubic will fit better than linear
prices = S0 - m.sqrt(S0 * K) * m.exp(-r*T)/np.pi * spline_cub( m.log(S0/K) )
return prices
# A class that stores option parameters (in order to write BS/Heston class neatly)
class Option_param():
def __init__(self, S0=10000, K=10000, T=.1, v0=0.04, payoff="call", exercise="European"):
"""
S0: current stock price
K: Strike price
T: time to maturity
v0: (optional) spot variance
exercise: European or American
"""
self.S0, self.v0, self.K, self.T, self.exercise, self.payoff = S0, v0, K, T, exercise, payoff
现在让我们使用 GEKKO：
#Initialize Model
m = GEKKO()
#define parameter
eq = m.Param(value=5)
#initialize variables
x1,x2,x3,x4,x5 = [m.Var(lb=-1, ub=1),m.Var(lb=1e-3, ub=1),m.Var(lb=1e-3, ub=1),m.Var(lb=1e-3, ub=20),m.Var(lb=1e-3, ub=1)]
#initial values
x1.value = 0
x2.value = 0.5
x3.value = 0.5
x4.value = 0.5
x5.value = 0.5
X = [x1,x2,x3,x4,x5]
#Equations,Feller Condition
m.Equation(2*x3*x4 - x2*x2 >=0)
def least_sq(x):
""" Objective function """
Heston_param = Heston_process(mu=0, rho=X[0], sigma=X[1], theta=X[2], kappa=X[3])
m = 1/(spreads**2)
if len(m) == 1:
l = 1
else:
l = (m - np.min(m))/(np.max(m)-np.min(m))
opt_param = Option_param(S0=s0, K=strikes, T=T, v0=X[4], exercise="European", payoff="call" )
Hest = Heston_pricer(opt_param, Heston_param)
prices_calib = Hest.FFT(strikes)
results = (l * (prices_calib-prices)**2)/len(prices)
return results
m.Obj(m.sum(least_sq(X)))
#Set global options
m.options.IMODE = 3 #steady state optimization
#Solve simulation
m.solve()
#Results
print('')
print('Results')
print('x1: ' + str(x1.value))
print('x2: ' + str(x2.value))
print('x3: ' + str(x3.value))
print('x4: ' + str(x4.value))
print('x5: ' + str(x5.value))
这里的问题是 ifft scipy 函数不起作用，因为 GEKKO.The 给出的变量类型不同问题是我不知道如何替换或避免它.
错误是这样的：
<ipython-input-16-305a2bb0769b> in fft_(K, S0, r, T, cf)
78
79 integrand = m.exp(- 1j * b *m.CV(np.arange(N))*dx) * cf(x - 0.5j) * 1/(x**2 + 0.25) * weight * dx/3
---> 80 integral_value = np.real(ifft(integrand)*N)
81 spline_cub = interp1d(ks, integral_value, kind="cubic") # cubic will fit better than linear
82 prices = S0 - m.sqrt(S0 * K) * m.exp(-r*T)/np.pi * spline_cub( m.log(S0/K) )
/usr/local/lib/python3.7/dist-packages/scipy/_lib/deprecation.py in call(*args, **kwargs)
18 warnings.warn(msg, category=DeprecationWarning,
19 stacklevel=stacklevel)
---> 20 return fun(*args, **kwargs)
21 call.__doc__ = msg
22 return call
<__array_function__ internals> in ifft(*args, **kwargs)
/usr/local/lib/python3.7/dist-packages/numpy/fft/_pocketfft.py in ifft(a, n, axis, norm)
274 a = asarray(a)
275 if n is None:
--> 276 n = a.shape[axis]
277 if norm is not None and _unitary(norm):
278 inv_norm = sqrt(max(n, 1))
IndexError: tuple index out of range
谁能帮我调试this.Thank你
Gekko 要求表达式不是黑盒，而是可以用特殊类型的变量（Gekko 类型）来表示，用于自动微分和稀疏检测。使用 Scipy.optimize.minimize 等求解器可能会更好地解决这个问题。这是一个comparison of the two on a simple problem.
Scipy
import numpy as np
from scipy.optimize import minimize
def objective(x):
return x[0]*x[3]*(x[0]+x[1]+x[2])+x[2]
def constraint1(x):
return x[0]*x[1]*x[2]*x[3]-25.0
def constraint2(x):
sum_eq = 40.0
for i in range(4):
sum_eq = sum_eq - x[i]**2
return sum_eq
# initial guesses
n = 4
x0 = np.zeros(n)
x0[0] = 1.0
x0[1] = 5.0
x0[2] = 5.0
x0[3] = 1.0
# show initial objective
print('Initial Objective: ' + str(objective(x0)))
# optimize
b = (1.0,5.0)
bnds = (b, b, b, b)
con1 = {'type': 'ineq', 'fun': constraint1}
con2 = {'type': 'eq', 'fun': constraint2}
cons = ([con1,con2])
solution = minimize(objective,x0,method='SLSQP',\
bounds=bnds,constraints=cons)
x = solution.x
# show final objective
print('Final Objective: ' + str(objective(x)))
# print solution
print('Solution')
print('x1 = ' + str(x[0]))
print('x2 = ' + str(x[1]))
print('x3 = ' + str(x[2]))
print('x4 = ' + str(x[3]))
盖柯
from gekko import GEKKO
m = GEKKO()
#initialize variables
x1,x2,x3,x4 = [m.Var(lb=1,ub=5) for i in range(4)]
x1.value = 1; x2.value = 5; x3.value = 5; x4.value = 1
m.Equation(x1*x2*x3*x4>=25)
m.Equation(x1**2+x2**2+x3**2+x4**2==40)
m.Minimize(x1*x4*(x1+x2+x3)+x3)
m.solve()
print('x1: ' + str(x1.value))
print('x2: ' + str(x2.value))
print('x3: ' + str(x3.value))
print('x4: ' + str(x4.value))