During the Cassini-Huygens mission, which lasted for nearly 20 years and came to an end by September 2017, the Saturnian moon Enceladus has been found to be one of the most interesting celestial bodies in our solar system regarding future exploration missions. Plumes that ascend from the south polar region could be detected by the Cassini spacecraft and were identified to contain ionized water vapor. Additionally, camera recordings showed large cracks in Enceladus’ kilometer thick icy shell with higher surface temperatures, the so-called “Tiger stripes”. An evaluation of the gathered information indicates cryovolcanism as the reason for such activity, leading to the conclusion that there likely is a global liquid water ocean below the surface . The possibility of proving extraterrestrial life or habitability was a starting signal for several subsequent mission concepts.
One of these ideas has been developed during the Enceladus Explorer (“EnEx”) collaborative project. A possible scenario would be to retrieve water samples from the ocean by drilling through the icy shell. If aiming for cracks close to the geysers and hitting liquid water reservoirs, the mission objective could be accomplished relatively easily. Therefore, a maneuverable melting probe has been built, the IceMole  (Fig. 1).
During several field tests, its ability to locate itself and to maneuver through ice has been shown. A scientific highlight was the first successful retrieving of a clean brine sample of the Blood Falls in Antarctica . After demonstrating the technical feasibility of the mission concept in principle, necessary technologies shall now be further enhanced in six follow-up projects, each focusing on different parts. In EnEx-CAUSE, a fully autonomous navigation system is developed, involving sensor fusion, decision making and steering of the probe.
Guiding the IceMole raises a sequence of challenging tasks. From modeling the system dynamics, to trajectory planning, up to determining the actuator controls, a coherent concept is required. A proven strategy to handle the problems of trajectory planning and model predictive control is given by the definition and solution of optimal control problems (OCP). Nowadays, complex and highly nonlinear systems do not necessarily represent an obstacle against the finding of feasible and even optimal solutions under real-time conditions. The combination of a binary character and the amount of 32 control variables, however, makes the IceMole a challenging candidate for real-time control. As a consequence, a methodology to deal with potentially high-dimensional mixed-integer control problems (MIOCP), that have to be solved within time restrictions, has to be devised.
I-a Related Work
The study of OCPs became popular in the 1950s and 1960s, when there was a high demand for solutions in the military sector . Solution strategies can be divided into two classes of approaches : On the one hand, direct methods transcribe the problem by collocation, multiple shooting or full discretization methods into nonlinear programs (NLP). These potentially high-dimensional problems can be solved by sequential quadratic programming (SQP) or interior point algorithms. On the other hand, the indirect approach uses Pontryagin’s maximum principle to state the necessary optimality conditions, which lead to a boundary value problem that again can be solved by shooting methods. MIOCPs, often also referred to as hybrid optimal control problems, have firstly been studied in the 1980s. They are considered to be much harder to solve than pure OCPs . For solving them, there is no universal best approach but a variety of strategies. By analogy with OCPs, they can be discretized to mixed-integer nonlinear programs (MINLP), which typically are solved by branch-and-bound algorithms. Alternatively, there are strategies such as relaxation, transformation or rounding techniques and combinations thereof, as e.g. proposed in  and .
Our contribution is twofold. Firstly, we present a methodology on how the IceMole melting probe can be autonomously steered during glacier tests using a combination of off-the-shelf optimization software. We propose a pragmatic concept that at some point neglects dynamics exactness, but guarantees feasible solutions under real-time conditions. By relaxing a MIOCP, we are able to plan and adjust optimal trajectories in a robust and efficient way. The binary variables are determined in a subsequent fitting algorithm. We set up a MINLP that is being frequently solved so that feasible and effective binary control configurations can be applied.
Secondly, we conduct a post-test analysis. Because of the limited repeatability of the test procedure, we do a comprehensive evaluation of the obtained data in order to enhance our algorithms for future field tests. On the one hand, this includes revealing and adjustment of weak aspects coming with our methodology. On the other hand, we perform automated parameter identification to improve our system dynamics model.
Ii System Dynamics
Ii-a Melting Probes
Melting probes are designed to reach a certain destination in an icy environment in order to perform a particular action. Typically, they allow only straight-downward motion. This type of movement has been studied in  and  with respect to the correlation between induced heating power and melting velocity under the effects of conductive loss. In order to actively maneuver through ice, a concept around curvilinear motion is required. For instance, there are studies about the use of heating plates. One  deals with the application of an asymmetrically distributed force, another  presents a framework around differential heating and a constant contact force. By applying a temperature gradient at the surface of a pressurized plate, a curvilinear melting path can be enforced, which can be described by the curve radius and the mean velocity. However, transferring these principles to a universal model for trajectory planning of a specific melting probe in three dimensions is not obvious. The inertia of the body cannot be disregarded, as well as specific attributes that come along with the IceMole’s construction.
Ii-B IceMole Model and Constraints
The IceMole is a special melting probe, as it uses two complementary advance strategies. A screw guarantees a steady contact force towards the ice front and, in principle, allows movement against the weight, while a collection of heater elements determines the actual speed and, additionally, allows curvilinear melting due to the differential heating concept.
In Fig. 2, the IceMole’s layout, including the positioning of the screw and all heater elements , is illustrated. The heater units are numbered clockwise, beginning with those disposed in the head at the front side , followed by those positioned in the walls and completed by those in the back part . Before we model the dynamic behavior in an efficient way, we make a simplification. The head heater units in each quadrant are grouped together, so that there are four groups . We attach one control variable to each of these groups. The same can be done with the wall heater units and the corresponding back heater unit on each side. Since the units in the back corners stick out of the body, they constantly have to be turned on (). Otherwise, the IceMole could get stuck or, at least, be restricted in its movement. Therefore, they are not considered as control variables for modeling purposes. Altogether, we define eight control variables .
To model the system dynamics, we do not want to examine the melting process itself, involving concepts of thermal conductivity, convection and the transfer of energy by phase changes. Instead, we will set up a relation between the observable state of the system and the actuator controls by using a priori knowledge of the system. Assuming that the IceMole’s translation can only occur along its longitudinal axis and neglecting gravitation, we describe our motion model by the kinematic equations of motion with six degrees of freedom in an east-north-up frame as:
It is stated as an autonomous ordinary differential equation system
. The state vectoris a pose, which is defined by coordinates and a quaternion . Quaternions allow a non-unique but differentiable representation of rotations in three dimensions.
It is still unclear how the heater controls are related to the kinematic equations. There is a general understanding of how each unit affects the system’s behavior. Namely, the IceMole’s speed is mainly defined by the front heaters, while its curvilinear motion is additionally determined by the wall heaters. The heaters on the back panel support backwards motion when the probe has to be pulled out of a partially frozen melting channel. During the EnEx collaborative project, a proven control strategy was to either operate on all head heaters for straight forward movement or to turn on half of the head heaters and the corresponding wall heaters for curvilinear melting. Since a nearly linear dependence between the melting velocity and the induced heat has been pointed out in , we model this relationship by
where is the maximal velocity.
Due to the torque applied by the screw, the IceMole holds a certain roll rate. Its extent depends on the ice quality, on the one hand, and on whether wall heaters are used or not, on the other hand. As previous tests have shown, those effects can change locally and cannot be characterized from data validation at this point. Therefore, we make a very simple generalization of a linear increase with respect to the velocity:
Assuming that we know the radius for maximal curvilinear melting and that rotation can only occur during forward movement, we define the angular velocities with respect to the y- and z-axis by
The parameter describes to which extent the curvilinear melting process can be attributed to the head and wall heater units respectively. In total, our model includes four parameters .
Power supply is restricted so that not all units can be operated at once. Therefore, power consumption for each group has to be incorporated, meaning that the weighted sum of the control variables has to be restricted:
where is the sum of the single units related to each group and is the maximal total power consumption minus the power, which has to be reserved for the back heater corner units . Additionally, the IceMole’s software interface demands that power has to be separately reserved for the wall and back heaters. This leads to one further constraint:
Iii Trajectory Planning
With the given dynamics, trajectories can be planned. Motion planning algorithms like PRM (Probabilistic roadmap) planners  or RRTs (Rapidly-exploring random trees)  are considered to be state of the art to formulate and solve these problems in many applications. Since we do not take any possible obstacles into account, the main difficulty consists in complying with the dynamics. The system of differential equations (1) has to be solved, as well as the constraint equations (2) - (7) for every time step. Additionally, one might want to plan trajectories with respect to different criteria. This favors a formulation as an optimization problem, allowing to quickly solve all at once.
Until this point, we have not looked closely at the possible control variables. The heater elements that are plugged into the IceMole can either be switched on or off, which means they should be modeled as binary controls. Under this assumption, the planning problem has to be stated as a MIOCP. Compared to OCPs, this lifts the planning problem to another level. Continuous optimal control problems can be solved for millions of variables and constraints in less than a second, while the computation time for MINLPs in principal grows exponentially with every additional optimization variable . Thus, we relax the MIOCP. The continuous heater controls will further be referred to as , while the binary controls are stated as . An OCP is defined in its standard form by
To solve a problem of this form, we use discretization methods, e.g. the trapezoidal rule to transcribe it into a NLP of general form
Common techniques to solve these sparse and potentially high-dimensional NLPs are sequential quadratic programming or interior point algorithms. We use the optimization routine WORHP  with an SQP method that uses an interior point algorithm on the QP level. Below, we will deduce our trajectory planning problem. We distinguish between two formulations in which one version is more restrictive and supposed to be solved only once, while the other one is executed frequently.
Iii-a Reference Trajectory
Reaching a fixed goal pose with an exact orientation is nearly impossible because of the IceMole’s limited controllability. Additionally, the final roll angle is of no particular interest. Therefore, only pitch and yaw angle shall be considered. The attitude deviation in the final state can be interpreted as a rotation. We use the boxminus operator () to determine its representation in Euler angles:
where are the deviations with respect to the three Euler angles. Now, can be restricted to a certain value :
Additionally, we introduce box constraints for the translational velocity and the angular velocities:
We define the values as control variables, since an increase in the number of linear equations is often numerically preferable over fewer nonlinear equations. Still unknown is the objective of the OCP. Since a linear dependence between the mean melting velocity and power consumption has been modeled, a compelling energy-optimal solution can hardly be obtained with respect to this model. Therefore, we choose a time-optimal formulation and add a weighted penalty term (with ) for the final orientation deviation. Altogether, with a state vector and a control vector , we define our trajectory planning problem by
To follow a computed path, a replanning or feedback control algorithm has to be implemented. If obstacles are considered, trajectory updates can involve large discrepancies. Since we neglect the occurrence of obstacles, the potential for error only arises from the dynamics, control application and environmental conditions. Therefore, we can expect a new solution to lie close enough to the previous one, which leads to a high convergence rate. However, the problem can become locally infeasible. Since the melting probe’s possibilities to change direction are rather limited, reaching a fixed destination might become more and more difficult with decreasing distance. Therefore, we remove the final position constraint and, instead, add another penalty term (with ) to the cost function. One could still bound this term, but we leave this decision to a higher-level decision making module. Given a current pose update , the trajectory replanning problem at time step can be stated as
It can be interpreted as a model predictive control algorithm with an infinite prediction and control horizon and weights not equal to zero only for the final state (cf. ).
Iv Actuator Control
Given the trajectories generated by MPC, a way to transform the real-valued control functions into binary control variables has to be found. This problem is closely related to pulse width modulation (PWM) techniques but with the restriction that signals cannot be separated and are connected through additional constraints. Sager  proposes sum up rounding strategies for either singular, independent control functions or those who are related by SOS1 constraints (). In our case, however, there are three inequalities of a general linear form, meaning that those strategies will not be sufficient. Therefore, we propose a MINLP optimization approach, which generates feasible control variables during every trajectory calculation cycle.
To determine the full set of binary controls , we ungroup the previously defined control variables:
Although we do not consider the dynamic equations at this point, we still have to account for their physical reasoning. If we look at a short time interval of adequate length , we demand that the defect for every control function , respectively , vanishes:
From now on, we inspect a replanning step of fixed length . The control function
is evaluated by applying a linear interpolation to the discretized solution ofMPC and is extended constantly over the integral bounds.
If we postulate term (13) to be minimal over the whole length for all control variables and introduce weights , this leads to the following cost function:
By applying an integration scheme, e.g. the chained trapezoidal rule, we obtain a quadratic objective. The binary control switching points are given by the supporting points for integration. Those share the same grid for the inner and the outer integral. If constant terms are neglected, the optimization problem can be stated as
We call it the binary control transformation (BCT). It is solved by using the MINLP software SCIP  with WORHP on NLP level. If the computation time is too long, which is not surprising in connection with MINLPs, the algorithm can be terminated and it still provides a feasible solution. If the replanning step size is small enough, the continuous control variables will merely change from one point to the other. Therefore, we can take the feasible solution and use it as an initial guess for the next calculation cycle:
As a result, we can expect solutions that become much better after a few iterations.
In August 2017, the EnEx initiative performed a field test on the Langenferner glacier in Italy, during which the IceMole was controlled fully automatically for the first time. The main objective was to show the general capability of all partners’ subsystems and proper communication between them. In two “long distance straight on” melting tests for around 25 meters each, the correct actuating of all heater units was proven. In terms of modeling dynamics and trajectory planning, however, those tests were less informative. The only way to verify this system’s capability to its full extent is given by curvilinear melting, which may be accompanied by serious complications when tested on a glacier. During the first EnEx project, there were only a few melting tests of this type in which the experience was gained that the IceMole could get stuck when being pulled backwards out of a curved melting channel. Therefore, we tested the maneuverability of the probe only for a relatively short-distance melting process of approximately three meters.
Before evaluation, we have to state to which terms the results could be obtained. The previously presented algorithms were performed in a sequence, as illustrated in Fig. 3. Two things have not been discussed before. At first, there is a software safeguard mechanism (SSG) that might prevent the IceMole from switching certain heater units on, e.g. due to overheating. Secondly, we did not have any ground truth information with respect to the position data. Data from different subsystems, i.e. the IMU (inertial measurement unit), the EnEx-RANGE APU (Autonomous Pinger Units) network  and the screw rotation counter, was gathered and continuously fused to a pose by the sensor fusion subsystem .
To analyze the position data, we make the following definition. From now on we consider trajectories from to as a sequence of states
In general, we can choose a representation defined by a parameter setting and a sequence of control variable sets . We can calculate a trajectory by forward integration by
For comparison, trajectories are generated by applying the integration scheme via the implicit trapezoidal rule to the different control set sequences, which have been computed and applied during the test run (cf. Fig. 3). The parameter selection was fixed to values , , and .
As can be deduced from Fig. 1, the IceMole has not been able to follow its planned trajectory very well. All trajectories that have been simulated predict the IceMole to move in a much narrower curve than in reality. We explain these distinctions by two major reasons:
approximation error introduced by the binary control transformation
inaccuracy of the system dynamics model
Below, we will locate these causes of error in detail and, consequently, improve our algorithms for future field tests.
Vi Post-Event Analysis
Vi-a Binary Control Transformation
With the BCT algorithm, a real-time capable strategy to determine and alternate 32 binary controls has been postulated. Since the algorithm does not necessarily terminate in time, optimality cannot be expected. More importantly, the system dynamics are neglected.
As can be deduced from Fig. 5, there is a systematic deviation regarding power distribution. During the test, we used weights for all controls within the objective function (14). A physically more reasonable choice would be to set . This would lead to an optimization of the heating power distribution. Besides, since we are solving a fitting problem with additional constraints, the method tends to converge to a solution from below. To enforce a higher power output to the detriment of the uniformity of the distribution, we add a second term with a comparatively small weight to our cost function (14):
Given these modifications, we can define the modified binary control transformation by
Additionally, there was an inconsistency regarding the power restrictions (6) and (7) between MPC and BCT, which has been eliminated. We performed a reoptimization with the modified algorithms and got a much better result, as can be seen in Fig. 6 with respect to the power distribution, and in Fig. 7 with respect to the trajectory.
|Power Consumption (rel. to MPC)|
Vi-B Parameter Identification
To a certain extent we have to accept that our model is not precise enough. We neither take any heating flow concepts into account nor kinetic equations. Additionally, we do not consider locally changing environmental conditions, e.g. small crevasses. Overall, we cannot expect physical correctness for a highly complex melting process if a model is that simple. Nevertheless, it should be of higher precision for this short-distance test case.
To overcome the inaccuracies, we will not extend our model in the first place, as its simplicity allows complex trajectory planning. Instead, we will adapt our system parameters. To recalculate a better choice, we state another optimization problem. With the forward integration scheme proposed in (17), a trajectory as a result of a known sequence of applied binary controls can be calculated with respect to parameters . To fit the fusion position data, we define a cost function by the least squares error with respect to the y- and z-position:
Finally, the parameter identification can be stated as an unconstrained optimization problem
We solve PI and obtain an optimal parameter set with , , and . Except for , this parameter set is reasonable with respect to its physical background. The comparatively high value of (around +20) can partially be explained by the effective loss of induced heating power along the process chain.
By modifying the BCT algorithm and using parameter identification, we resolved the different sources of error separately. However, on the basis of the given data set, we cannot bring these parts together without deliberating. The challenge would be to find a parameter set that is feasible at planning time. Because of the step-wise decrease of induced heating power along the process sequence, the previously computed parameter set is not directly applicable. To overcome this, there are two possible approaches. One approach would be to integrate BCT and SSG into the MPC algorithm. As to BCT, it is too computationally expensive and contradicts the separation of a MIOCP in the first place. SSG is not an accessible algorithm. Therefore, both would have to be approximated with respect to their effects. Alternatively, one could try to “invert” the applied controls at first, in order to perform PI for a continuous control set . Again, this requires an approximation of the inverse operator. Both strategies could be carried out by introducing parameter-dependent functions that model the assumed decrease of induced heating power along the sequence of algorithms.
Vii Conclusions & Future Work
We presented an optimization-based framework for the maneuvering of an autonomous melting probe. It includes a strategy for reformulating and separating a MIOCP into two consecutive algorithms, so that an altogether strong solution can be found under real-time conditions. The framework’s fundamental applicability was proven during an intermediate glacier test. An evaluation of the data indicated systematic errors which could be attributed to different sources. Those defects were fixed separately by using parameter identification and adapting an objective function. A combination of these strategies yields a new problem which will become a subject to further investigations.
As a result, we have to adapt our framework so that the cycle between trajectory (re-)planning and parameter identification can be closed. We will improve our models with the objective to integrate an on-line adaptive model for the next glacier test.
This work was supported by the German Aerospace Center (DLR) with financial means of the German Federal Ministry for Economic Affairs and Energy (BMWi), project “EnEx-CAUSE” (grant No. 50 NA 1505).
-  C. Porco, P. Helfenstein, P. Thomas, A. Ingersoll, J. Wisdom, R. West, G. Neukum, T. Denk, R. Wagner, T. Roatsch et al., “Cassini Observes the Active South Pole of Enceladus,” Science, vol. 311, no. 5766, pp. 1393–1401, 2006.
-  B. Dachwald, J. Mikucki, S. Tulaczyk, I. Digel, C. Espe, M. Feldmann, G. Francke, J. Kowalski, and C. Xu, “IceMole: a maneuverable probe for clean in situ analysis and sampling of subsurface ice and subglacial aquatic ecosystems,” Annals of Glaciology, vol. 55, no. 65, pp. 14–22, 2014.
-  J. Kowalski, P. Linder, S. Zierke, B. von Wulfen, J. Clemens, K. Konstantinidis, G. Ameres, R. Hoffmann, J. Mikucki, S. Tulaczyk et al., “Navigation Technology for Exploration of Glacier Ice With Maneuverable Melting Probes,” Cold Regions Science and Technology, vol. 123, pp. 53–70, 2016.
-  H. J. Sussmann and J. C. Willems, “300 Years of Optimal Control: From the Brachystochrone to the Maximum Principle,” IEEE Control Systems Magazine, vol. 17, no. 3, pp. 32–44, 1997.
-  O. Von Stryk and R. Bulirsch, “Direct and indirect methods for trajectory optimization,” Annals of Operations Research, vol. 37, no. 1, pp. 357–373, 1992.
-  C. Kirches, Fast Numerical Methods for Mixed-Integer Nonlinear Model-Predictive Control. Springer, 2011.
-  M. Gerdts, “A variable time transformation method for mixed-integer optimal control problems,” Optimal Control Applications and Methods, vol. 27, no. 3, pp. 169–182, 2006.
-  S. Sager, Numerical methods for mixed-integer optimal control problems. Der andere Verlag Tönning, Lübeck, Marburg, 2005.
-  H. Aamot, Heat Transfer and Performance Analysis of a Thermal Probe for Glaciers, ser. Technical report. U.S. Army Materiel Command, Cold Regions Research & Engineering Laboratory, 1967.
-  S. Ulamec, J. Biele, O. Funke, and M. Engelhardt, “Access to glacial and subglacial environments in the solar system by melting probe technology,” Reviews in Environmental Science and Bio/Technology, vol. 6, no. 1-3, pp. 71–94, 2007.
-  H. Kumano, A. Saito, S. Okawa, and Y. Yamada, “Direct contact melting with asymmetric load,” International Journal of Heat and Mass Transfer, vol. 48, no. 15, pp. 3221–3230, 2005.
-  K. Schüller, J. Kowalski, and P. Råback, “Curvilinear melting–a preliminary experimental and numerical study,” International Journal of Heat and Mass Transfer, vol. 92, pp. 884–892, 2016.
-  R. F. Stengel, Flight dynamics. Princeton University Press, 2015.
-  L. E. Kavraki, P. Svestka, J.-C. Latombe, and M. H. Overmars, “Probabilistic roadmaps for path planning in high-dimensional configuration spaces,” IEEE transactions on Robotics and Automation, vol. 12, no. 4, pp. 566–580, 1996.
-  S. M. Lavalle, “Rapidly-exploring random trees: A new tool for path planning,” Tech. Rep., 1998.
-  R. E. Bixby, M. Fenelon, Z. Gu, E. Rothberg, and R. Wunderling, “Mixed integer programming: a progress report,” The Sharpest Cut, pp. 309–325, 2004.
-  C. Büskens and D. Wassel, “The ESA NLP Solver Worhp,” in Modeling and Optimization in Space Engineering. Springer, 2012, pp. 85–110.
-  M. Bloesch, H. Sommer, T. Laidlow, M. Burri, G. Nuetzi, P. Fankhauser, D. Bellicoso, C. Gehring, S. Leutenegger, M. Hutter et al., “A primer on the differential calculus of 3d orientations,” arXiv preprint arXiv:1606.05285, 2016.
-  L. Grüne and J. Pannek, “Nonlinear model predictive control,” in Nonlinear Model Predictive Control. Springer, 2011, pp. 43–66.
-  S. J. Maher, T. Fischer, T. Gally, G. Gamrath, A. Gleixner, R. L. Gottwald, G. Hendel, T. Koch, M. E. Lübbecke, M. Miltenberger, B. Müller, M. E. Pfetsch, C. Puchert, D. Rehfeldt, S. Schenker, R. Schwarz, F. Serrano, Y. Shinano, D. Weninger, J. T. Witt, and J. Witzig, “The SCIP Optimization Suite 4.0,” ZIB, Takustr.7, 14195 Berlin, Tech. Rep. 17-12, 2017.
-  D. Eliseev, D. Heinen, K. Helbing, R. Hoffmann, U. Naumann, F. Scholz, C. Wiebusch, and S. Zierke, “Acoustic in-ice positioning in the enceladus explorer project,” Annals of Glaciology, vol. 55, no. 68, pp. 253–259, 2014.
-  J. Clemens, “Multi-robot in-ice localization using graph optimization,” in Autonomous Robot Systems and Competitions (ICARSC), 2017 IEEE International Conference on. IEEE, 2017, pp. 36–42.