我在做天气数据模拟，有四个方程，包含四个变量。有很多其他参数，但它们的值是提供的。我正在使用gekko optimizer并得到这个错误。我是编程新手。请帮忙。我已尽力提供尽可能多的信息。

# all libraries

from scipy.integrate import quad
import numpy as np

# all input parameters are as below

h = 1000 # height of air slab in meters
ps = 100000 # surface presuure in pascals
ph = 88000 # pressure at slab top in pascals
atop = 0.2 # entrainment parameter
b = 0.3 # moistening parameter
zh = 0.2 # hydrologically active soil depth in meters
zt = 0.4 #  thermally active soil depth in m
n = 0.25 # soil porosity
vmin = 0.5 # Mineral volume fraction of soil
vorg = 0.25 # Organic volume fraction of soil
c = 1 # Exponent in evaporation efficiency
eta = 1 # Coefficient in runoff ratio
r = 2 # Exponent in runoff ratio
uz = 4 # Mixed layer wind speed in m/s
l = 500000 # Length scale of region
qin = 8/1000 #Effective humidity of incoming air
g = 9.81 # gravitational accelaration in m/s**2
rd = 287.058 # gas constant for dry air in SI uits
cp = 1005 # dry air specific heat at constant pressure in SI units
sigma = 5.670374419*10**(-8)
A = 0.75 # integration term
m = 1/7 # integration term
c1 = 0.001 # coefficient for forced velocity
c2 = 0.00025 # coefficient for buoyancy velocity
LHV = 22.5*10**5 # latent heat of vaporization of water in SI units
rv = 461 # gas consatant for water vapor in SI units
rho = 1.225 # air density in SI units
RS = 1000 # incoming solar radiation in SI units
yh = (ph/ps)**(3/2)

# calculation of incoming water vapor flux

Qin = ((ps-ph)/g)/l*uz*qin
def f1(y):
    return (y**(8/3*rd/cp))*((y-yh)**(m-1))
int1 = quad(f1, yh, 1)
def f2(y):
    return (y**(8/3*rd/cp))*((1-y)**(m-1))
int2 = quad(f2, yh, 1)
from gekko import GEKKO
m = GEKKO()
s = m.Var(value = 0.5)
qm = m.Var(value = 1)
tg = m.Var(value = 250)
thetam = m.Var(value = 350)
m.Equation(((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))-(((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))*eta*(s**r))-((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))==0)
m.Equation(Qin-(((ps-ph)/g)/l*uz*qm)+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))-((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))==0)
m.Equation(RS*(1-(0.17-0.085*s))+(1.24*((qm*ph)/(0.622*thetam*(ph/ps)**(rd/cp)))**(1/7)*sigma*(thetam*(ph/ps)**(rd/cp))**(4))(1-(0.75*((2/3*qm*ps/g*(1-(ph/ps)**(3/2)))**(1/7))))+(sigma*(thetam**(4))*m*A*((2/3*qm*ps/g)**(m))*int2)-(sigma*(tg**(4)))-((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*(tg-thetam)*rho*cp)-LHV*((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))==0)
m.Equation(((1.24*((qm*ph)/(0.622*thetam*(ph/ps)**(rd/cp)))**(1/7)*sigma*(thetam*(ph/ps)**(rd/cp))**(4))+(sigma*(tg**(4))))*0.75*((2/3*qm*ps/g*(1-(ph/ps)**(3/2)))**(1/7))-sigma*(thetam**(4))*m*A*((2/3*qm*ps/g)**(m))*int2-sigma*thetam**(4)*m*A*(2/3*qm*ps/g)**(m)*int1+1.2*(c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*(tg-thetam)*rho*cp==0)
m.solve()