Model predictive control for electrodynamic tether geometric profile in orbital maneuvering with finite element state estimator

This paper studies the control of the geometric profile of a librating electrodynamic tether by model predictive control using the induced electric current in tether only. First, a high-fidelity multiphysics model of an electrodynamic tether system is built based on the nodal position finite element method and the orbital-motion-limited theory. Second, a state estimator is proposed to estimate the geometric profile of a librating electrodynamic tether, where only the positions and velocities at the tether ends are measurable. The non-measurable geometric profile of the tether between two ends is estimated by the high-fidelity multiphysics model with the input of the measurement at tether ends in the spatial domain. The extended Kalman filter is applied to estimate the geometric profile of the tether while avoiding the singularity or ambiguity in the estimation. Third, the geometric profile control of a librating electrodynamic tether is converted into a trajectory tracking problem of the underactuated electrodynamic tether system. The induced electric current in the tether is the only control input. Its value is optimized by the model predictive control method subject to the output and input control constraints. The numerical simulation results show that the proposed approach can effectively control the shape of the librating electrodynamic tether to the reference trajectory.


Introduction
Space tethers, especially the bare electrodynamic tethers (EDT), have been proposed for removing space debris from Earth orbits for decades because of low mass, compact size, propellant-free, and ease of operation [1][2][3][4]. However, the application of the EDT is impeded by its intrinsically unstable libration motion [2], which is caused by the continuous interaction of the current-carrying EDT with the Earth's magnetic field [3,5]. Many efforts have been devoted to stabilizing the libration motions of an EDT by regulating the electric current in tether as the only control input [5][6][7][8]. The tether was first modeled by a simple dumbbell model to apply analytical control laws [8,9] with a constant electric current assumption [10]. Later, a flexible tether model was proposed to discretize the tether into a finite number of segments with an energy-rate feedback control law [7]. Here, the geometric profile of the tether is assumed measurable by the Kalman filter. Unfortunately, it is practically impossible to place sensors along the thin tether to measure its spatial position. Nevertheless, the discretized tether model results in multiple sets of (inplane and out-of-plane) libration angles instead of only one set of libration angles in the dumbbell model, see Fig. 1 where the dashed lines represent the dumbbell model. This renders the existing libration control strategies based on the dumbbell model not applicable. To mitigate this problem, the authors proposed a virtual straight tether that connects the first and last notes of the flexible tether model to represent the libration motion, see Fig. 1, and developed a libration stabilization strategy by an electric current on-off control scheme to keep the Hamiltonian of the virtual tether bounded [5,11]. Although effective, this approach works only if the real tether is straight or close to straight. Once the EDT is bent noticeably by orbital perturbative forces, the Hamiltonian of the virtual tether no longer represents the real EDT, and the energy-based libration control scheme fails [11].
Motivated by the challenge, this work develops a novel approach based on the finite element method, which can estimate the geometric profile of tether by the input of positions and velocities at the tether ends only and then to keep it straight. First, the spatially continuous dynamical system is discretized into finite elements (FE). The FE model establishes the relationship in the spatial domain between the measured state at the tether ends and the non-measurable state along the tether. Then, the extended Kalman filter is applied to propagate the full state (measured and non-measurable) in the time domain to avoid the singularity if the degrees of freedom of the non-measurable state are greater than that of the measured state. Random noises are introduced in the measurement model to account for the sensor noises, unmodeled dynamics, and model uncertainties associated with the finite element method [12]. Once the entire state estimation is obtained, the model predictive control (MPC) method is applied to control the geometric profile of EDT by minimizing the difference between the real and virtual straight tethers, as shown in Fig. 1. The finite element discretization is completed by the nodal position finite element method (NPFEM) [3,5], where the states are the nodal positions and the orbital height and tether length are calculated from the nodal positions. Thus, the possible round-off errors in the computation caused by the huge different scale between the orbital height and tether length are effectively avoided. Finally, the symbols used in this paper is given in Table 1. 2 High-fidelity dynamic model of electrodynamic tether

System description
Consider an EDT orbiting the Earth as shown in Fig. 2. It consists of a main and a sub-spacecraft Fig. 1 Schematic diagram of a discretized EDT with virtual straight tethers connected by an elastic tether that is bare and electrically conductive [5]. Electric potential bias is induced in the tether as the latter interacts with the Earth's magnetic field. Electrons in ambient plasma are attracted by a positively charged segment of the bare EDT [2,13] and emitted back to the plasma by an emitter at the cathodic end [13], resulting in an electric current in the EDT. The current-carrying EDT interacts with the Earth's magnetic field and results in a Lorentz (electrodynamic) force against EDT's orbit motion.
Four orthogonal coordinate systems are defined to describe the motion of an EDT [11,14], as shown in Fig. 2. The first is the global inertial coordinate system (oxyz) with its origin o at the Earth center, x-axis in the vernal equinox, z-axis aligned with the Earth's rotation axis, and y-axis forming a right-hand coordinate system. The second is the Earth-fixed coordinate system (ox f y f z f ) where the gravity field, atmosphere, and ionospheric plasmas densities of the Earth are described [14]. The ox f y f z f is assumed to coincide with the oxyz initially. The third is the orbital at the center of mass (CM) of the EDT system, x o -axis from the center of the Earth to the CM, y o -axis in the orbital motion direction, and z o -axis forming a right-hand system. Finally, the local coordinate system (o l tnb) is defined at the element level. The origin o l is at the CM of each element, and the detailed definition can be found in [11,14,15].

High-fidelity dynamic model of EDT
A high-fidelity multiphysics model of the EDT has been built, and the tether is spatially discretized by the nodal position finite element method [16,17]. To make the paper self-contented for interested readers, this model is briefly described here. The main and subspacecraft are simplified as lumped masses due to the large ratio of the tether length over the size of the spacecraft. In the current study, the tether is discretized into (n-1) two-noded elastic tensile elements with n nodes. The dynamic model of the discretized EDT is derived from the principle of virtual work [11,17], where the overdot denotes the time derivative, x t ð Þ2 R 3n is the vector of nodal position, (m and k t ð Þ2 R 3nÂ3n ) are the mass and stiffness matrices, and c t ð Þ ¼ a 1 m þ a 2 k t ð Þ ½ 2R 3nÂ3n is the Rayleigh damping with a 1 and a 2 being damping coefficients [12,16]. It is noted the mass matrix m is a constant matrix, (a) (b) Fig. 2 The schematic diagram of a flexible EDT model which is advantageous in numerical integration [16,17]. The lumped masses of spacecraft (m m and m s ) are added directly to the diagonal elements mass matrix m corresponding to the first and last nodes. Their values are given in the simulation section. The terms (f g t ð Þ, f d t ð Þ, and f l t ð Þ2 R 3n ) are the vectors of gravity, atmosphere drag, and electrodynamic force, respectively.
The gravity and atmosphere drag are evaluated by, where the subscript j 2 ½1; n À 1 denotes the j th element, S is the matrix of linear shape function [17], g g t ð Þ ¼ T f 2g t ð Þg f t ð Þ is the vector of gravity acceleration evaluated by the Earth's gravitational model [11] with superscripts g and f representing the global inertial and Earth-fixed coordinate systems, f l d is the vector of atmospheric drag per unit length with the superscript l representing the local frame, T f 2g t ð Þ and T l2f t ð Þ represent the transformation matrix between different frames: ox f y f z f 7 !oxyz and o l tnb7 !ox f y f z f [5,11], respectively, q a is evaluated by the Naval Research Laboratory Mass Spectrometer Incoherent Scatter Radar (NRLMSIS-00) model [14], (c d ¼ 2:2 for the tether, c d ¼ 2:0 for the spacecraft), _ x l and v l a t ð Þ ¼ v l;t a ; v l;n a ; v l;b a À Á T are the velocity vectors of the EDT and velocity of surrounding atmosphere with superscript l representing the local coordinate system [17]. The electrodynamic force vector is evaluated by, where I j t ð Þ is the induced electric current of the j th element. It is assumed constant within element, is the Lorentz force per unit length with superscript l representing the local frame o l tnb, b l j t ð Þ is the vector of magnetic field strength evaluated at the element's CM by the International Geomagnetic Reference Field model (2014) [11,14], t is a unit vector along the element's taxis in the local coordinate frame [11,14], d t ð Þ 2 R 3nÂm is the mapping matrix from the induced electric current u t ð Þ 2 R m to the force vector f l t ð Þ2 R 3n . As shown in Fig. 2b, the bare EDT can be divided into a positively biased segment AB that collects electrons and a negatively biased segment BC that collects ions [4,13], the induced electric current I s; t ð Þ and potential bias U s; t ð Þ obey the orbital-motionlimited (OML) theory [4,13,18], Þare the coefficients relate to the electron density of ambient plasma and tether physical properties [4,18], and e m s; t ð Þ is the electric motive force (EMF) [14].
The boundary conditions that reflect the design of the EDT system are chosen for the maximum efficiency of the current generation for deorbit purpose [18], such that, where l B is the length of the positively biased segment of the tether and is unknown in priori [18], l C is the total length of EDT, [4,18].
Equations (10)-(12) form a nonlinear two-point boundary value problem. They are discretized by the FE method by the same mesh as in Eq. (1) [4,18]. If the point B is coincided with the pp th node (pp = 2 * n-1), see Fig. 2b, the discretized form of the OML theory is written as where subscripts pos and neg denote the positive and negative potential bias properties, h is the coefficient matrix, r edt is the vector of nodal electric current and potential bias, and f emf pos and f emf neg are the EMF vectors with subscripts ''pos'' and ''neg'' representing the positively and negatively charged segments of the tether [18].
Correspondingly, the boundary condition equation (12) is rewritten as, where is inside an element, its location is first approximated by the solution of OML theory with the boundary conditions at two ends of the element. Then, the element is divided into two sub-elements and modeled by Eqs. (13)- (14). Detailed information can be found in our previous work [3,18].
Finally, the discretized OML theory of an EDT system is assembled by a standard FE procedure, where h edt 2 R 2nÂ2n , r edt ¼ I 1 ; U 1 ; . . .I pp ; U pp ; . . .; À I n ; U n Þ T 2 R 2n and f emf 2 R 2n [18]. Equation (17) is a nonlinear algebraic equation with unknown boundaries at two ends and must be solved by the Newton-Raphson algorithm iteratively [3,4,18] to obtain the electric current along the tether at each node, see where I EDT 2 R n . Recast Eq. (1) in the state-space form for MPC development, where the superscript c denotes the continuous time system, E 1 2 R 3nÂ3n is an identity matrix, and the control input (electric current) u t ð Þ 2 R m must satisfy the following inequality, where u low t ð Þ and u upp t ð Þ2 R m are the lower and upper bounds of control input (the induced electric current). For the deorbit application, u low t ð Þ = 0 and u upp t ð Þ ¼ GI EDT with G 2 R mÂn given in (8). As demonstrated in [3,18], the value of u upp t ð Þ depends on the libration motion, the geometric profile, physical properties of the EDT tether, and environmental parameters. Thus, u upp t ð Þ changes dynamically at each time step and is updated by solving the discretized OML equation at the beginning of each time step.
3 Model predictive control for geometric profile of electrodynamic tether

State estimator
The entire state of the tether is divided into two parts: the measured state at two ends of the EDT and the nonmeasurable state inside the tether. In the following, we will derive a state estimator like a virtual sensor estimating the non-measurable state inside the tether, based on the input of measured state at the ends (boundary) of the tether in the temporal and spatial domains. First, linearize the kinematics of the EDT system in the state space as the system model, where the superscript d denotes the discrete-time system, the subscripts obs and (k-1 and k) denote the system model at the times k and k-1, respectively, W obs;kÀ1 is the noise model accounting for the model and parameter uncertainties, N 0; Q obs À Á represents the normal distribution with the zero mean 0 and covariance matrix Q obs 2 R 6nÂ6n , and u kÀ1 2 R m is the vector of input control at time k-1.
Then, the state measurement model is constructed in two parts: a measurement model of the measurable state located at the boundary of the tether and a virtual measurement model of the non-measurable state inside the tether.
For the measurable state Y 1 obs 2 R q 1 ¼12 at the boundary of the tether qX, the measurement model is defined as where the superscript 1 denotes the measurement model of measurable state, H 1 obs 2 R q 1 Â6n is the measurement matrix, and V 1 obs;kÀ1 2 R q 1 is the measurement noise of the sensors, such that, Here, E 2 2 R 3Â3 is an identity matrix.
For the non-measurable state X 2 obs 2 R 3nÀq 1 2 X inside the tether, the virtual measurement model is derived from the finite element model of the EDT [11,17] in the spatial domain with the input of measured state Y 1 obs 2 R q 1 ¼12 at the boundary qX, such that, where the superscript 2 denotes the measurement model of non-measurable state, Y 2 obs;kÀ1 2 R 6nÀq 1 is the estimated state, h 2 X obs;kÀ1 À Á 2 R 6nÀq 1 is a nonlinear virtual measurement function based on the spatially discretized finite element model of the flexible EDT, and V 2 obs;kÀ1 2 R 6nÀq 1 is the noise model that includes the noise propagated from the input of measured state at the boundary, the discretizing error of the finite element method and unmodeled dynamics of the EDT, G kÀ1 is the constraint equation of the first measurement model as listed in Eq. (27), and oG kÀ1 =ox obs is the Jacobian matrix of the first measurement model with respect to the state x obs .
The V 1 obs;kÀ1 and V 2 obs;kÀ1 are assumed to be independent zero-mean Gaussian white-noise process, such that, where d i;j and EX Á; Á ½ are the Kronecker delta and expectation functions, respectively, the subscripts i and j represents the i and j component of a vector, R 1 obs 2 R q 1 Âq 1 and R 2 obs 2 R 6nÀq 1 ð Þ Â6nÀq 1 ð Þ are the known noise covariance matrices. Equation (30) is numerically solved by the Lagrangian multiplier method directly without reducing the degree of freedom of the system [12]. However, the solution Y 2 obs is not unique if 6n À q 1 ð Þ[ q 1 , i.e., the dimension of the non-measurable state vector Y 2 obs 2 X is greater than that of the measurable state vector Y 1 obs at the boundary oX. To avoid the ambiguity or singularity in the reconstruction of the non-measurable internal state, the extended Kalman filter is used to propagate the full state of the EDT X k with the known initial condition and the hybrid measurement model This state estimator is defined as the finite element Kalman filter (FEKF) estimator that contains the prediction and correction stages [19][20][21] as follows.
• Prediction step: • Correction step: where the superscripts -or ? denote prior or posterior estimators, K obs is the Kalman gain matrix, is the Jacobian matrix of the nonlinear measurement model with H 2 k ¼ oh 2 =oX obs j X obs ¼X À obs;k being the Jacobian matrix of the nonlinear virtual measurement function h 2 ðÞ, H 1 is a constant matrix as listed in Eq. (28), R obs ¼ diag R 1 obs ; R 2 obs À Á , and E 3 2 R 6nÂ6n is an identity matrix.
The stability of the newly proposed FEKF estimator can be proved based on the direct method of Lyapunov for the extended Kalman filter. The detailed proof can be found in [19,20]. The observability of the newly proposed FEKF estimator is proved by the Popov-Belevitch-Hautus (PBH) test, see the page 174 of [22], and is briefly outlined here. The pair A c obs ; H in the continuous system equation (22) and the measurement model is observable only and only if when the PBH matrix is a full column rank for all eigenvalues k eigen 2 ÉR, such that, where A c obs ¼ diag E 1 ; 0 ð Þis the state transition matrix with the superscript c representing the continuous time system.
Since the term (k eigen E 1 À A c obs ) becomes singular only at the eigenvalues of the observed system, it suffices to merely investigate the observability of the system modes [22,23]. The measurement matrix Á is known, and only the diagonal elements exist. Obviously, the rank of the PBH matrix of the state-space model remains as the full column matrix.

Model predictive control
Once the state estimator is established, the MPC method will be employed to optimize the input control sequence considering constraints, such as state, output, and control (the induced electric current) for each element of a flexible EDT model explicitly [24][25][26][27].
Consider the high-fidelity model of EDT in the state-space form as listed in (22). It is discretized in the time domain by the zero-order hold method, such that, where X k ¼ X þ obs;k is the state of the tether estimated by the state estimator in Sect. 3 k with e ðÞ and Dt being the exponential function and sample time interval, respectively. In the current work, the exponential function e ðÞ is calculated by the Dgpadm algorithm [28], and the inverse of a matrix ðÞ À1 is calculated by the pseudoinverse algorithm.
By performing the recurrence computation of (39)-(40), the predicted state sequence X kþ1 can be written in a compact form as a function of the state X k and the input state sequence u k at the current time k . .
The values of n P and n c will be given at the simulation section.
Substituting (41) into (42) leads to The control objective is to minimize the difference between the actual state, or geometrical profile, of the tether Y t ð Þ 2 R 3n and the desired or reference state Y ref t ð Þ 2 R 3n of a virtual straight tether that connects the main and sub-spacecraft, subject to the system equation in (47) and the constraints Here, the cost function J k at the current time k is constructed based on the standard MPC procedure [29,30], u k t ð Þ is the controlled input, Q mpc 2 R 3nÂ3n and P mpc 2 R mÂm are the output error weight and input control weight matrices, respectively, u low where , and S T mpc Q mpc S mpc is a constant term that does not affect the minimization of the cost function [25,31]. The superscript T represents the transpose of a matrix.
Rewriting the input constraint in (51) in terms of the vector u k leads to The state constraint in (52) is introduced to prevent a large departure of the predicted state Y kþ1 from the previous state to ensure a smooth control performance. This is achieved by defining the lower and upper bounds of the predicted vector as Y low kþ1 ¼ 0:95Y kþ1 and Y upp kþ1 ¼ 1:05Y kþ1 . Based on Eq. (47), Eq. (52) is reformulated in terms of the control vector u k as, where the matrix D 2 2 R 2 n c Âm ð Þ and the vector P 2 2 R 2 n c Âm ð Þ are formed as follows, Combining Eqs. (54) and (56) by concatenating one on top of the other leads to a single constraint equation, where

Constraint softening
Considering the fact that the EDT system is underactuated because there is only one control input (the electric current along the tether), the constraint in Eq. (58) may occasionally be breached. Instead of imposing the hard restriction of the constraint, which may lead to the infeasible hard-constrained problem, a constraint-softening technique is employed to penalize the constraint violation. Based on the slack variable method [32,33], the problem in Eqs. (53) and (58) can be rewritten as, where e 1 2 R q 2 is the vector of slack variables denoting the constraint violation, P 4 2 R q 2 Âq 2 is the positive-definite weight matrix for the slack variable vector e 1 , and e upp 1 is a vector of the upper bound of the slack variables. When the base MPC problem is infeasible, the partial slack variables in Eq. (59) become positive to allow the constraints relaxed, which makes the problem feasible.
Introduce a new vector u 0 k ¼ u k ; e 1 ð Þ T . Substituting the new vector into Eq. (59) yields, with where E 5 2 R q 2 Âq 2 is an identity matrix.

Generation of reference trajectory for the tether state
The desired or reference trajectory of the tether state is a straight line that passes through the main and subspacecraft. To ensure the length of the reference state is the same as the real tether, a special procedure is needed to map the nodes along the real tether to the reference tether see Fig. 3. It contains two steps as follows, (i) Generate the end nodes of reference tether First, generate a straight line from the main spacecraft (point A) to pass the sub-spacecraft (point C) and end at point C 0 . The total length of the straight reference tether AC 0 shall be the same as the instant length of the real tether l t ð Þ, such that, (ii) Generate the internal nodes along the reference tether After the reference trajectories of position are obtained, the reference trajectories of internal nodes of the real tether are mapped to the reference tether based on the instant length of each element one by one starting from the main spacecraft. Then, the reference trajectories of all nodes of the reference tether are obtained, By using the state space equation (39), the reference trajectory Y kþ1 ¼ Y kþ1 k j ; Y kþ2 k j ; À Á Á Á ; Y kþn c k j ; Á Á Á ; Y kþn p k j Þ T can be easily obtained.

Stability analysis of proposed MPC problem
The asymptotic stability of MPC is examined through the monotonicity condition of cost function, The monotonicity condition of the cost function of a constrained MPC is guaranteed if any of the following terminal conditions is applied (i) terminal equality constraints, (ii) terminal set constraints, and (iii) terminal penalty function [31,32]. It is known that the terminal penalty function is superiority to the first two methods due to its wider feasible region [31,34]. Accordingly, recasting the constrained MPC in Eq. (60) by adding a terminal inequality constraint, such that, where e 2 2 R 3n and Q N p 2 R 3n is the penalty weight matrix that should satisfy the monotonicity condition of cost function [32], such that where K np is the optimal gain of a linear quadratic regulator controller with the same Q mpc and P mpc Directly Measurable node Unmeasurable node Fig. 3 Schematic diagram for generation of reference tether weighting matrices, which is obtained by solving the corresponding Riccati equation [31,32].

Implementation of proposed MPC with FEKF estimator
The proposed MPC with the FEKF estimator for controlling the geometric profile of a flexible EDT in the deorbit process is shown in Fig. 4. The implementation process is shown in Fig. 5. The constrained MPC problem is translated into a quadratic programming (QP) problem and solved by a QP solver.

Numerical simulation
The model is built up and simulated at the global inertial frame as listed in Eq. (1). The results are then transformed from the global frame to the orbital frame to make it consistent with the research work in space engineering [11,14]. The formulation of the orbital frame can be found in our previous works, where the detailed conversion process can be found.

System parameters and initial conditions
The effectiveness of the proposed MPC algorithm with the FEKF estimator is demonstrated by numerical simulation. It is applied to control the geometric profile of an EDT in the deorbit process by regulating the induced electric current I s; t ð Þ in tether with two initial tether profiles, as shown in Fig. 6, where the tether is initially straight (case A) and bent into a halfcircle (case B). As a comparison, each case is analyzed with and without control input by the MPC method. The physical properties of the EDT are shown in Table 2. The electrical boundary conditions, which represent the EDT system design, are given in Table 3. The standard deviation R 1 obs and R 2 obs for the measurement model in (32) is assumed as 20 m as per [35]. The tether is uniformly divided into five elements, and the time step is set as 0.001 s. The libration angles of the reference state are used to represent the libration motion of the discretized EDT model.
The prediction and control horizons are selected as n p = 20 and n c = 5, and the weight matrices are set as P mpc ¼ 1000E 2 and Q mpc ¼ 0:1E 2 to make the tracking error as a primary optimization goal, based on a priori knowledge [25,36] and the trade-off between tracking error and control effect.
The lower and upper bounds of the input current in Eq. (51) are defined as zero and the maximum possible induced electric current, respectively. This choice of the upper bound of input current is to maximize the electron collection efficiency of the bare tether [4,18]. First, the dynamic simulation is conducted without any control by (1), where the electric current is obtained by numerically solving (17) [18]. The numerical simulation stops once the libration angles of the EDT exceed 90 degrees, which implies the tether flips over. The simulation results are shown in Figs. 8, 9, and 10. Second, the dynamic simulation is conducted with the MPC in (60). The simulation results are plotted in Figs. 7, 8, and 9 for comparison. Figure 7 shows the in-plane libration angle of the EDT with/without MPC quickly exceeds 90°( &1.57 rad) because the MPC does not include the constraints on the libration angles in the current work. The EDT flips over and becomes unstable see Fig. 8. This is expected because there is no control for the libration angle. The Lorentz force keeps pumping energy into the libration motion of the EDT system. However, it is interesting to see in Fig. 8 that the profile of the EDT is much closer to the straight virtual tether when the proposed MPC is applied. The application of MPC results in a longer stable time before flipping over than the case without electric current control, see Fig. 8. The electric current at the cathodic end of EDT is shown in Fig. 9. It shows the induced electric current in the EDT is at the maximum level in the case without control as designed, see Fig. 9a. The profile of the electric current with MPC is shown in Fig. 9b, where the profile of the electric current is clearly regulated between the minimum value (zero) and the upper bound. This explains why the profile of tether with control is closer to the straight tether, and the stable period is about 15% longer than the case without electric current control.

Case B
In this case, the tether is initially bent into a half-circle shape, and both ends are in the local vertical. Similarly, two simulations are conducted with and without MPC until the tether flips over. Simulation results are shown in Figs. 10, 11, and 12. Figure 10 shows the libration motion of the EDT for Case B. It clearly shows that the in-plane angle of EDT without tether profile control quickly flips over 90°as expected due to no libration angle control. On the contrary, the in-plane angle of the EDT with tether profile control is significantly less when the proposed controller is applied. For instance, the in-plane angle of EDT with the electric current regulation by the MPC is less than 0.2 rad, while the in-plane angle of EDT without the electric current control reaches 90°in the same period. This demonstrates that the proposed MPC effectively straightens the tether toward the straight reference state in the local vertical initially and results in a much smaller in-plane libration angle, see Fig. 10b and 11. Finally, Fig. 12 shows the electric current profile at the cathodic end is significantly different in cases with and without MPC because the corresponding tether profiles are quite different. This is different from Case A, where the profiles of tether in both controlled and uncontrolled conditions are Table 2 Physical properties of EDT [4] Name  similar. In summary, the proposed MPC with the FEKF estimator can straighten the flexible bent tether in the libration motion.

Conclusion
In this paper, a model predictive control (MPC) with finite element Kalman filter (FEKF) estimator is proposed to straighten a bent flexible electrodynamic tether or prevent a straight electrodynamic tether from bending in the deorbit process. This is achieved by controlling the induced electric current in the tether. The non-measurable tether profile is effectively obtained by the proposed FEKF estimator with the input of the measured state at both ends of the tether. By tracking a virtual straight tether that passes through both ends of the tether as the reference state, the framework of the finite element is based on the model predictive control algorithm, which includes the orbital-motion-limited theory and constraints on the input electric current and the control errors, is developed. Moreover, both the stability of the proposed MPC and the observability of the state estimator are examined analytically. The problem of controlling the geometric profile of the electrodynamic tether is transferred into a trajectory tracking problem and then solved by nonlinear programming. Numerical case studies successfully demonstrate that the proposed model predictive control method with FEKF estimator effectively controls the geometric profile of the electrodynamic tether toward the straight reference state with a single control input of electric current. Finally, it is noted there are two limitations of the proposed method, which could be improved in future. First, the maximum libration angles of the EDT can be introduced to the constraints of the state to limit the libration angles. Second, the desired/reference trajectory of the virtual straight tether varies at each time step, which make it hard to track under the current constraint of the control input. In the future, the proposed control strategy can be compared with the simple current switch on-off control in the deorbit process by the EDT system.