Providing dynamic efficiency threshold as function

Idea

We can give the efficiency threshold as a function with the signature:

def eta_r(t,oe,mass,battery):
    # ... compute eta as a function of the inputs ...
    return eta_r

For example, this can be a function based on the battery level:

def eta_r(t,oe,mass,battery):
    return min(400/(4*battery - battery_dod), 1)

Demonstration

We first do our imports and define the initial/final state and weights

import numpy as np
from numpy.random import rand
import matplotlib.pyplot as plt
import time

import sys
sys.path.append("../../../")    # path to pyqlaw
import pyqlaw

# initial and final elements: [a,e,i,RAAN,omega,ta]
LU = 42164.0
GM_EARTH = 398600.44
VU = np.sqrt(GM_EARTH/LU)
TU = LU/VU

rp_gto = 200 + 6378
ra_gto = 35786 + 6378
sma_gto = (rp_gto + ra_gto)/(2*LU)
ecc_gto = (ra_gto - rp_gto)/(ra_gto + rp_gto)
KEP0 = [sma_gto,ecc_gto,np.deg2rad(23),0,0,0]
KEPF = [1,0,np.deg2rad(3),0,0,0]
oe0 = pyqlaw.kep2mee_with_a(np.array(KEP0))
oeT = pyqlaw.kep2mee_with_a(np.array(KEPF))
woe = [3.0, 1.0, 1.0, 1.0, 1.0]

# spacecraft parameters
MU = 1500
tmax_si = 1.0    # 1 N
isp_si  = 1500   # seconds
mdot_si = tmax_si/(isp_si*9.81)  # kg/s

# non-dimensional quantities
mass0 = 1.0
tmax = tmax_si * (1/MU)*(TU**2/(1e3*LU))
mdot = np.abs(mdot_si) *(TU/MU)
tf_max = 10000.0
t_step = np.deg2rad(15)

Suppose we want to define eta as a function of battery; we define some battery related quantities

# battery levels
battery_initial = 3000*3600/TU            # Wh --> Ws --> W.TU
battery_dod = 500*3600/TU
battery_capacity = (battery_dod,battery_initial)
charge_rate = 1500          # W
discharge_rate = 500        # W
battery_charge_discharge_rate = (charge_rate, discharge_rate)
require_full_recharge = True

We now construct the QLaw object, and set the problem with information about the battery.

# construct problem
prob = pyqlaw.QLaw(
    integrator="rk4",
    elements_type="mee_with_a",
    verbosity=2,
    print_frequency=500,
    use_sundman = True,
)

# set problem
prob.set_problem(oe0, oeT, mass0, tmax, mdot, tf_max, t_step,
    battery_initial = battery_initial,
    battery_capacity = battery_capacity,
    battery_charge_discharge_rate = battery_charge_discharge_rate,
    require_full_recharge = require_full_recharge,
    woe = woe)
prob.pretty()

Transfer:
  a  : 5.7800e-01 -> 1.0000e+00 (weight: 3.00)
  f  : 7.3009e-01 -> 0.0000e+00 (weight: 1.00)
  g  : 0.0000e+00 -> 0.0000e+00 (weight: 1.00)
  h  : 2.0345e-01 -> 2.6186e-02 (weight: 1.00)
  k  : 0.0000e+00 -> 0.0000e+00 (weight: 1.00)

We now define the functions eta_a() and eta_r() that will be called at each step of the integration. Here, we make use of the relative efficiency only, so the eta_a() function is merely a placeholder (we might as well pass eta_a=0.0 later), but is still defined as a function for the sake of demonstration.

# efficiency functions
def eta_a(t,oe,mass,battery):
    return 0.0    # not using absolute efficiency threshold

def eta_r(t,oe,mass,battery):
    return min(400/(4*battery - battery_dod), 1)

Finally, we call the solve() method of prob, passing the functions eta_a and eta_r as arguments

# solve
tstart_solve = time.time()
prob.solve(eta_a=eta_a, eta_r=eta_r)
prob.pretty_results()

 iter   |  time      |  del1       |  del2       |  del3       |  del4       |  del5       |  el6        |
      0 |  3.201e-02 | -4.2200e-01 |  7.3009e-01 |  0.0000e+00 |  1.7727e-01 |  0.0000e+00 |  6.4336e-01 |
    500 |  5.813e+01 | -3.8563e-01 |  6.7137e-01 |  1.4516e-03 |  1.3048e-01 |  3.1424e-04 |  1.2788e+02 |
   1000 |  1.224e+02 | -3.3722e-01 |  6.0152e-01 |  5.2874e-03 |  8.4814e-02 |  8.0152e-04 |  2.5382e+02 |
   1500 |  1.949e+02 | -2.7639e-01 |  5.1111e-01 |  1.1664e-03 |  4.3824e-02 | -8.6243e-04 |  3.7922e+02 |
   2000 |  2.819e+02 | -1.3341e-01 |  3.7391e-01 |  1.8435e-03 |  1.6334e-02 | -2.5591e-04 |  5.0430e+02 |
   2500 |  3.976e+02 |  1.0253e-02 |  1.1446e-01 |  4.9404e-03 |  2.3246e-03 | -7.1226e-05 |  6.2619e+02 |
Target elements successfully reached!
Exit code : 2
Converge  : True
Final state:
  a  : 1.0013e+00 (error: 1.2623e-03)
  f  : 2.5225e-03 (error: 2.5225e-03)
  g  : 7.8086e-03 (error: 7.8086e-03)
  h  : 2.6280e-02 (error: 9.3977e-05)
  k  : 1.0635e-04 (error: 1.0635e-04)
Transfer time : 448.94083725747794
Final mass    : 0.8640773711998154

and we can visualize the results

# plots
fig1, ax1 = prob.plot_elements_history(to_keplerian=True, TU=TU/86400, time_unit_name="day")
fig2, ax2 = prob.plot_trajectory_3d(sphere_radius=0.1)
fig3, ax3 = prob.plot_controls()
fig4, ax4 = prob.plot_battery_history(TU=TU/86400, BU=TU/3600,
    time_unit_name="day", battery_unit_name="Wh")
fig5, ax5 = prob.plot_efficiency(TU=TU/86400, time_unit_name="day")
plt.show()

Using 4489 steps for evaluation