CODEWITHSUNDEEP
pythonoptimizationcontrolsgekko
trimat
trimat

Reputation: 101

How to limit the SoC in Optimal control?

I'm trying to build this physiological model. I want to model athlete recovery, this is similar to an Electric Car state of charge. I'm unsure how to limit the 100% SoC.

The issue I'm having is that the problem gets very large since I need to get the previous SoC to measure the current Soc. Furthermore, I'm not sure how to model the fact that the discharge rate is not the same as the charge rate. Thank you

#discharging
W'Exp = (P - CP) * dt                     
soc = soc - W'Exp                       
#If the Power for the data point is below CP or CV, W' is reconstituted:
W'Rec = (CP - P) * dt                            
Rate = (W' - W'Bal) / W'                         
soc = soc + W'Rec * Rate

Upvotes: 1

Views: 116

Answers (2)

John Hedengren
John Hedengren

Reputation: 14366

Check out a Gekko implementation of the simple energy storage model found here as an example problem. You may be able to modify this model for the athlete.

periodic

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as mtick
from gekko import GEKKO

m = GEKKO(remote=False)

t = np.linspace(0, 24, 24*3+1)
m.time = t

m.options.SOLVER   = 1
m.options.IMODE    = 6
m.options.NODES    = 3
m.options.CV_TYPE  = 1
m.options.MAX_ITER = 300

p = m.FV()           # production
p.STATUS = 1
s = m.Var(100, lb=0) # storage inventory
store = m.SV()       # store energy rate
vy = m.SV(lb=0)      # store slack variable
recover = m.SV()     # recover energy rate
vx = m.SV(lb=0)      # recover slack variable

eps = 0.7

d = m.MV(-20*np.sin(np.pi*t/12)+100)

m.periodic(s)

m.Equations([p + recover/eps - store >= d,
             p - d == vx - vy,
             store == p - d + vy,
             recover == d - p + vx,
             s.dt() == store - recover/eps,
             store * recover <= 0])
m.Minimize(p)

m.solve(disp=True)

#%% Visualize results
fig, axes = plt.subplots(4, 1, sharex=True)

ax = axes[0]
ax.plot(t, store, 'C3-', label='Store Rate')
ax.plot(t, recover, 'C0-.', label='Recover Rate')

ax = axes[1]
ax.plot(t, d, 'k-', label='Electricity Demand')
ax.plot(t, p, 'C3--', label='Power Production')

ax = axes[2]
ax.plot(t, s, 'C2-', label='Energy Inventory')

ax = axes[3]
ax.plot(t, vx, 'C2-', label='$S_1$')
ax.plot(t, vy, 'C3--', label='$S_2$')
ax.set_xlabel('Time (hr)')

for ax in axes:
    ax.legend(bbox_to_anchor=(1.01, 0.5), \
              loc='center left', frameon=False)
    ax.grid()
    ax.set_xlim(0, 24)
    loc = mtick.MultipleLocator(base=6)
    ax.xaxis.set_major_locator(loc)

plt.tight_layout()
plt.show()

Upvotes: 0

Jay
Jay

Reputation: 61

For the first problem, I guess SOC is a variable. The upper and lower bound of a variable can be set as follows:

m = GEKKO()    
SOC = m.Var(10, lb=0, ub=100)

10 is the initial value of SOC. For your second problem, are P, CP, W'Bal given?

Upvotes: 1

Related Questions