Search code examples
python-3.xdynamicmodelingnonlinear-optimizationgekko

Defining control variables and objective function in GEKKO


I've the following dynamical system that I am trying to solve in GEKKO

(1) d ϕ d t = − M T D M ϕ

I followed examples given here to set up the following equations in GEKKO

M, MT = get_mmt()
A = MT @ np.diag(Dhat) @ M
A = A[1:ngrid - 1]

# first node
m.Equation(phi_hat[0].dt() == 0)

# interior nodes
int_value = -A @ phi_hat  # function value at interior nodes
m.Equations(phi_hat[i].dt() == int_value[i] for i in range(0, ngrid-2))

# terminal node
m.Equation(phi_hat[ngrid-1].dt() == Dhat[end] * 2 * (phi_hat[end-1] - phi_hat[end]))

The state variables are phi_hat. I want to minimize the squared error difference between phi_meas and phi_hat from the model. Therefore, the following setting are used.

phi_hat = [m.Var(value=phi_0[i]) for i in range(ngrid)]
phi_hat = m.CV(value=phi_meas)
phi_hat.FSTATUS = 1  # fit to measurement phi obtained from 'def actual'

The error has to be minimized by optimizing parameter D (same as Dhat below) in equation 1.

Dhat0 = 500*np.ones(ngrid-1)
Dhat = m.FV(value=Dhat0[i] for i in range(0, ngrid-1))
Dhat.STATUS = 1  # adjustable parameter

I hope these settings are right. Please correct me if these are wrong. I am a little confused here, not sure if I should be using Dhat = m.CV(value=Dhat0[i] for i in range(0, ngrid-1)) in place of m.FV. Dhat is a vector of control variables that have to be estimated.

Post this, I'd like to ask for help on how to define the objective function (f = sum((phi(:) - phi_tilde(:)).^2)) using m.Minimize().

I have got the following settings for simulation

m.options.IMODE = 5  # simultaneous dynamic estimation
m.options.NODES = 5  # collocation nodes
m.solve()

EDIT: Detailed description of the problem that I am trying to solve can be found here. Note: Equation 1 mentioned here is actually equation 2 here.


Solution

  • For your question about using FV or CV:

    From your description, it sounds like it should be an FV. As a first step in solving your problem, I recommend that you verify that you have the specifications correctly set by switching to m.options.IMODE=7 and solve. This will tell you if you have the same number of equations and variables and also initializes the model so that parameter estimation is faster.

    If you use a CV then the squared error objective function is defined for you with m.options.EV_TYPE=2. You can also use m.Minimize() if you define the CV as a Var instead.