Train length and etcs

mseetc/ocp.py

@@ -6,7 +6,7 @@
 
 from mseetc.track import computeDiscretizationPoints
 
-from mseetc.utils import Options, var, postProcessDataFrame, splitLosses
+from mseetc.utils import Options, var, postProcessDataFrame, splitLosses, computeTunnelFactor
 
 
 class OptionsCasadiSolver(Options):
@@ -27,6 +27,10 @@ def __init__(self, paramsDict):
 
         self.integrateLosses = False  # integrate losses or take mid-point rule
 
+        self.chooseClosestTunnelCrossSection = True # if exact tunnel cross section is not specified in train tunnel resistances, choose the closest value
+
+        self.pwcLengthDependentTrackAttributes = False
+
         super().__init__(paramsDict)
 
 
@@ -121,7 +125,7 @@ def __init__(self, train, track, optsDict={}):
 
         # track parameters
 
-        self.points = computeDiscretizationPoints(track, numIntervals)
+        self.points = computeDiscretizationPoints(track, numIntervals, opts)
         self.steps = np.diff(self.points.index)
 
         # real-time parameters
@@ -193,16 +197,28 @@ def __init__(self, train, track, optsDict={}):
 
                 # gradient and curvature of current index
                 grad = self.points.iloc[i]['Gradient [permil]']/1e3
+                gradLinearTerm = self.points.iloc[i]["Gradient linear term [permil/m]"]/1e3
                 curv = self.points.iloc[i]['Curvature [1/m]']
+                curvLinearTerm = self.points.iloc[i]["Curvature linear term [1/m^2]"]
+                crossSection = self.points.iloc[i]['CrossSection [m^2]']
+                tunnelFactor = computeTunnelFactor(crossSection, train, opts)
 
                 # acceleration constraints
-                g += [trainModel.accelerationFun(ca.vertcat(time[i], velSq[i]), ca.vcat(u), grad, curv)]
+                g += [trainModel.accelerationFun(ca.vertcat(time[i], velSq[i], 0), ca.vcat(u), grad, gradLinearTerm, curv, curvLinearTerm, tunnelFactor)]
                 lbg += [accMin]
                 ubg += [accMax]
 
                 # coupling constraints
-                out = trainIntegrator.solve(time=time[i], velocitySquared=velSq[i], ds=self.steps[i],
-                    traction=Fel[i], pnBrake=Fpb[i], gradient=grad, curvature=curv)
+                out = trainIntegrator.solve(time=time[i],
+                                            velocitySquared=velSq[i],
+                                            ds=self.steps[i],
+                                            traction=Fel[i],
+                                            pnBrake=Fpb[i],
+                                            gradient=grad,
+                                            gradientLinearTerm=gradLinearTerm,
+                                            curvature=curv,
+                                            curvatureLinearTerm=curvLinearTerm,
+                                            tunnelFactor=tunnelFactor)
 
                 xNxt1 = ca.vertcat(time[i+1], velSq[i+1])
                 xNxt2 = ca.vertcat(out['time'], out['velSquared'])
@@ -322,7 +338,8 @@ def solve(self, terminalTime, initialTime=0, terminalVelocity=1, initialVelocity
         # initial guess
         # NOTE: good idea vel0 to be compatible with f0 (power-wise) to avoid nans at first iteration
 
-        vel0 = (60/3.6)**2
+        vel_avg = (self.points.index[-1]-self.points.index[0])/(terminalTime-initialTime) * 0.95
+        vel0 = vel_avg*vel_avg
         dt = (terminalTime - initialTime)/self.numIntervals
         t0 = initialTime
 
@@ -437,4 +454,4 @@ def solve(self, terminalTime, initialTime=0, terminalVelocity=1, initialVelocity
 
     else:
 
-        print("Solver failed!")
+        print("Solver failed!")