使用 m.time 随时间更新 gekko 变量及其先前值

Question

考虑一个涉及 MV、m.Var、m.CV 和 m.Param 矩阵的 gekko 方程，其中包括基于非线性回归的最大化问题。

如何在当前代码中添加一个与时间相关的决策变量“covar”，这样：

covar（在时间 t）= [1，gekko_attr1 ** attr1_pow，gekko_attr2 ** attr2_pow，gekko_attr3 ** attr3_pow] * covar（在时间 t-1）* 转置（[1，gekko_attr1 ** attr1_pow，gekko_attr2 * * attr2_pow, gekko_attr3 ** attr3_pow]) + 常量

数据集：

https://drive.google.com/file/d/10tVHghTvRThsDfyzprINdlVTtJTUbSyQ/view?usp=sharing

def gekko_fun():

    xm1 = np.array(data['Attr1'])
    xm2 = np.array(data['Attr2'])
    xm3 = np.array(data['Attr3'])
    ym = np.array(data['Result'])

    m = GEKKO(remote=False)

    # parameters
    gekko_attr1 = m.Param(value = xm1)
    gekko_attr2 = m.Param(value = xm2)
    gekko_attr3 = m.Param(value = xm3)

    intercept = m.MV(lb=0, name="Intercept")
    attr1_multiplier = m.MV(lb=0)
    attr2_multiplier = m.MV(lb=0)
    attr3_multiplier = m.MV(lb=0)

    intercept_lambda = m.MV(lb=0, ub = 1)
    attr1_lambda = m.MV(lb=0, ub = 1)
    attr2_lambda = m.MV(lb=0, ub = 1)
    attr3_lambda = m.MV(lb=0, ub = 1)

    attr1_pow = m.MV(lb=0, ub=1)
    attr2_pow = m.MV(lb=0, ub=1)
    attr3_pow = m.MV(lb=0, ub=1)
    
    intercept.STATUS = 1
    attr1_multiplier.STATUS = 1
    attr2_multiplier.STATUS = 1
    attr3_multiplier.STATUS = 1
    
    intercept_lambda.STATUS = 1
    attr1_lambda.STATUS = 1
    attr2_lambda.STATUS = 1
    attr3_lambda.STATUS = 1

    attr1_lambda_pow.STATUS = 1
    attr2_lambda_pow.STATUS = 1
    attr3_lambda_pow.STATUS = 1
    
    control_value = m.CV(value=ym)
    control_value.STATUS=1

    pred_value = m.Var
    error = m.Var

    maximized_value = m.Var

    m.options.IMODE = 2
    m.options.MAX_ITER = 999

    m.Equation(pred_value == (
        np.dot(
            [[1, gekko_attr1**attr1_pow, gekko_attr2**attr2_pow, gekko_attr3**attr3_pow]]
            ,#multiply
            np.add(
                np.dot(
                    [
                        [intercept_lambda,0,0,0],
                        [0,attr1_lambda,0,0],
                        [0,0,attr2_lambda,0],
                        [0,0,0,attr3_lambda]
                    ]                                                  
                    ,#multiply
                    [
                        [intercept],
                        [attr1_multiplier],
                        [attr2_multiplier],
                        [attr3_multiplier]
                    ]
                ),#add
                np.dot(
                    [
                        [attr1_multiplier,attr2_multiplier,attr3_multiplier],
                        [attr1_multiplier,0,0],
                        [0,attr2_multiplier,0],
                        [0,0,attr3_multiplier]
                    ],#multiply
                    [
                        [gekko_attr1**attr1_pow],
                        [gekko_attr2**attr2_pow],
                        [gekko_attr3**attr3_pow]
                    ]
                )
            )
        )
    ))
    
    m.Equation(error==control_value-pred_value)

    m.Equation(maximized_value == -m.log(covar) - covar*(error)**2)
    
    m.Maximize(maximized_value)

我是 gekko 的新手，如果我能在解决这个问题上得到一些帮助，我将不胜感激。 感谢您的宝贵时间。

Answer 1

如果存在离散时间方程而不是连续时间方程（例如微分方程），请使用 IMODE=3（默认）并使用数组构建变量。这是x[t] = x[t-1]+1 的示例问题。

from gekko import GEKKO
m = GEKKO()
x = m.Array(m.Var,10)
for i in range(9):
    m.Equation(x[i+1]==x[i]+1) 
m.Minimize((x[9]-10)**2)
m.solve()
print(x)

解决办法是：

 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :   5.399999994551763E-003 sec
 Objective      :   0.000000000000000E+000
 Successful solution
 ---------------------------------------------------
 
[[1.0] [2.0] [3.0] [4.0] [5.0] [6.0] [7.0] [8.0] [9.0] [10.0]]

Gekko 数组也可以进行 Numpy 矩阵运算。这是一个例子：

from gekko import GEKKO
import numpy as np
m = GEKKO(remote=False)
ni = 3; nj = 2; nk = 4
# solve AX=B
A = m.Array(m.Var,(ni,nj),lb=0)
X = m.Array(m.Var,(nj,nk),lb=0)
AX = np.dot(A,X)
B = m.Array(m.Var,(ni,nk),lb=0)
# equality constraints
m.Equations([AX[i,j]==B[i,j] for i in range(ni) \
                             for j in range(nk)])
m.Equation(5==m.sum([m.sum([A[i][j] for i in range(ni)]) \
                                    for j in range(nj)]))
m.Equation(2==m.sum([m.sum([X[i][j] for i in range(nj)]) \
                                    for j in range(nk)]))
# objective function
m.Minimize(m.sum([m.sum([B[i][j] for i in range(ni)]) \
                                 for j in range(nk)]))
m.solve()
print(A)
print(X)
print(B)

这是另一个例子：

import numpy as np
from gekko import GEKKO

m = GEKKO(remote=False)

# Random 3x3
A = np.random.rand(3,3)
# Random 3x1
b = np.random.rand(3,1)
# Gekko array 3x3
p = m.Array(m.Param,(3,3))
# Gekko array 3x1
y = m.Array(m.Var,(3,1))

# Dot product of A p
x = np.dot(A,p)
# Dot product of x y
w = np.dot(x,y)
# Dot product of p y
z = np.dot(p,y)
# Trace (sum of diag) of p
t = np.trace(p)

# solve Ax = b
s = m.axb(A,b)
m.solve()

对于未来的问题，请发帖Minimal, Reproducible example 以更好地诊断问题。