Effect of Reaction Wheel Imbalances on Attitude and Stabilization Accuracy

In this paper, the study of stabilization accuracy of a satellite equipped with a set of reaction wheels (RW) is presented. The model of motion takes into account possible static and dynamic reaction wheel imbalances. Due to the complexity of the model, which leads to the numerical issues, the effects of dynamic and static imbalances on inertial stabilization are studied analytically. As a result, estimations of the attitude and stabilization accuracy are presented in closed form.


Introduction
Modern small satellites are able to solve many scientific and applied problems, from conducting measurements of the Earth's magnetosphere to remote sensing. Most of these problems require precise satellite pointing and stabilization, which is usually provided by gyroscopic attitude control systems such as Reaction Wheels and Control Moment Gyros. These actuators, thanks to miniaturization, now could be installed even on-board the nanosatellites like 3-6U CubeSats [1][2][3][4].
RWs offer rather good performance and accuracy. However, there are several aspects that should be taken into account. The first is the saturation problem, which can be solved by auxiliary actuators such as magnetorquers. The second problem is the vibrations affecting spacecraft hardware, its flexible parts and overall dynamics [5,6], which might be crucial for some missions. One of the vibration sources is the RW imbalance. Typical angular rate of RW is several thousand rotations per minute, therefore, even the small imperfections in their balancing might significantly affect the performance of the attitude control system. There are two types of imbalances that are usually distinguished. The first is the static one, which appears when the RW center of mass is not located at its rotation axis. The second is the dynamic imbalance, which corresponds to the misalignment between the RW principal axis of inertia and the rotation axis (see Figure 1).

Introduction
Modern small satellites are able to solve many scientific and applied problems, from conducting measurements of the Earth's magnetosphere to remote sensing. Most of these problems require precise satellite pointing and stabilization, which is usually provided by gyroscopic attitude control systems such as Reaction Wheels and Control Moment Gyros. These actuators, thanks to miniaturization, now could be installed even on-board the nanosatellites like 3-6U CubeSats [1][2][3][4].
RWs offer rather good performance and accuracy. However, there are several aspects that should be taken into account. The first is the saturation problem, which can be solved by auxiliary actuators such as magnetorquers. The second problem is the vibrations affecting spacecraft hardware, its flexible parts and overall dynamics [5,6], which might be crucial for some missions. One of the vibration sources is the RW imbalance. Typical angular rate of RW is several thousand rotations per minute, therefore, even the small imperfections in their balancing might significantly affect the performance of the attitude control system. There are two types of imbalances that are usually distinguished. The first is the static one, which appears when the RW center of mass is not located at its rotation axis. The second is the dynamic imbalance, which corresponds to the misalignment between the RW principal axis of inertia and the rotation axis (see Figure 1).

Preliminary Remarks
Paper [25] describes the coupled model of motion. We present the brief derivation of a similar model here. It has a different final form that turned out to be more convenient for the software implementation and analytical study.
The spacecraft consists of the main hull with several RWs attached. The hull and each RW are supposed to be rigid bodies with known mass and inertia properties. Satellite position is described by point O-it is the arbitrary fixed point of the hull. Each reaction wheel is described by point O k (any point of the RW rotation axis) and axis of rotation e k (see Figure 2).

Preliminary Remarks
Paper [25] describes the coupled model of motion. We present the brief derivation of a similar model here. It has a different final form that turned out to be more convenient for the software implementation and analytical study.
The spacecraft consists of the main hull with several RWs attached. The hull and each RW are supposed to be rigid bodies with known mass and inertia properties. Satellite position is described by point O-it is the arbitrary fixed point of the hull. Each reaction wheel is described by point k O (any point of the RW rotation axis) and axis of rotation k e (see Figure 2).
The following reference frames are used in the paper. Inertial Frame is fixed in space (e.g., it might correspond to J2000 coordinate system). Body Frame is located at the hull point O, and its axes are fixed with respect to the hull, k-th RW Frame is attached to the RW in such a way that its third axis corresponds to the rotation axis. In order to derive the equations of motion the general equation of dynamics [28] for the system with ideal constraints is used: Here, the summation is for all points of the system (for rigid bodies, summation is replaced by integration), l is the system point index, l m is the point mass, l R  is its acceleration, l F is the total force affecting the point, l δ R is the point virtual displacement.
Every point of the hull i R and RW kj R is described by (see Figure 1):  The following reference frames are used in the paper. Inertial Frame is fixed in space (e.g., it might correspond to J2000 coordinate system). Body Frame is located at the hull point O, and its axes are fixed with respect to the hull, k-th RW Frame is attached to the RW in such a way that its third axis corresponds to the rotation axis.
In order to derive the equations of motion the general equation of dynamics [28] for the system with ideal constraints is used: Here, the summation is for all points of the system (for rigid bodies, summation is replaced by integration), l is the system point index, m l is the point mass, .. R l is its acceleration, F l is the total force affecting the point, δR l is the point virtual displacement.
Every point of the hull R i and RW R kj is described by (see Figure 1): where R O is the satellite radius-vector, r i is the vector from point O to the hull point, ρ k = OO k , ρ kj is the radius-vector from point O k to RW point. As was mentioned earlier, O k is the arbitrary point of RW rotation axis. It is reasonable to choose it as the projection of RW center of mass at its rotation axis. In this case, if the center of mass is on the rotation axis (no static imbalance) then O k represents the RW center of mass. Index i indicates the Here, δR O , δθ correspond to the hull virtual displacements, δϕ k is the k-th RW infinitesimal rotation. Displacements δR O , δθ, δϕ k are independent and correspond to N = 6 + n system degrees of freedom (n is the RW's quantity, and 6 for translational and rotational degrees of freedom for a rigid hull).
Accelerations of each satellite and RW points are: ..
where ω is the satellite angular velocity with respect to the Inertial Frame, Ω k = .
ϕ k e k is the RW angular velocity with respect to the hull. Finally, Equation (1) can be rewritten as follows: Here, M int k is the internal torque generated in the k-th RW axis that consists of the control and friction torques.

Equations of Motion Derivation
In order to derive equations of motion, it is necessary to gather all terms for the same virtual displacement and equate them to zero (since they are independent). This leads to the following equations: These equations can be simplified. The exact derivation is rather bulky and presented in the Appendix A. The final form of the motion equations is provided below.
Let us introduce the cross-product matrix: [a] × := The hull center of mass position with respect to the point O is given by: Similarly, the center of mass of k-th RW with respect to the point O k is: Total satellite mass m = m s + ∑ k m k . Denote the hull tensor of inertia with respect to the point O: and k-th RW tensor of inertia with respect to the point O k : Total tensor of inertia of the system with respect to the point O is: The total forces acting on the hull and RW are: Torques of the forces affecting RW with respect to the attachment point O k are: Similarly, the torque affecting the hull with respect to the point O is: To decrease the number of brackets in the equations, the following rule is used: a 1 × a 2 × . . . × a n ≡ a 1 × (a 2 × (. . . × a n ) . . .) Full satellite state vector is described by the hull position R O , velocity V O = .
R O , attitude quaternion Q = q 0 q T , Q = q 2 0 + q·q = 1 (optionally can be replaced by different attitude representation, e.g., Euler angles or direction cosine matrices), angular velocity ω, current RW rotation angle ϕ k and its angular velocity Ω k = . ϕ k . Satellite kinematics are: Equations of the satellite dynamics are: Here S is the following symmetrical matrix: where I k = e T k I k e k is the RW moment of inertia along its rotation axis, E 3×3 is 3 × 3 identity matrix, and N is: Note that the hull is considered to be a rigid body, so J s and r s remain constant in the Body Frame. However, in general case, ρ kc , I k and J depend on the current RW rotation angle ϕ k , and will change in the Body Frame. In addition, satellite acceleration (i.e., . V O ) is usually calculated in the Inertial Frame, while attitude dynamics are calculated in the Body Frame, which must be taken into account during the simulation.

Important Special Cases
The first special case to be considered is when there is no so-called static imbalance of RW, i.e., its center of mass is located at the rotation axis. In this case, it is reasonable to choose O k as RW center of mass. Therefore, O under Constraint (2) would denote the center of mass of the whole system, and ρ kc = 0. System Kinematics (3) remains the same, and dynamics becomes: If, in addition, there is no dynamical imbalance (i.e., the rotation axis is the RW dynamical symmetry axis), then I k (ω × Ω k ) + Ω k × I k ω = 0, I k e k = I k e k , e T k (ω × I k ω) = 0. This leads to the following equations of motion: In the matrix form: The simplified conventional model of the satellite motion is (see e.g., [29]): It is rather similar to the simplest Model (6), with the major difference in the equation that describes the RW rotation: it is actually affected by satellite rotation and external torques, e.g., by the gravity gradient torque. Note that internal torques from friction and control devices are usually much larger than external ones, so the effect of these additional terms is negligible.

Verification of the Motion Model
The model presented in Section 2 is rather complex. The conservation laws are used to validate its software implementation. In case of no external forces and torques, the total momentum in the Inertial Frame, angular momentum with respect to the system center of mass in Inertial Frame, and kinetic energy of the system must preserve, i.e., they are first integrals.
Total momentum of the system is: Velocities of the satellite points are: Hence, using Constraint (2), we obtain: Total angular momentum with respect to the system center of mass is: where R C is the system center of mass. Under Constraint (2) it can be written as follows: In the case when RW center of mass is at the rotation axis R C = R O . Simplification yields: Total kinetic energy is: After simplification: The numerical simulation is carried out to demonstrate the behavior of the abovementioned integrals with different integration steps. The system parameters are presented in Table 1. RW parameters are taken from [30]. The set of three identical RWs is considered. Their nominal axes of symmetry are aligned with the hull principal axes of inertia and each center of mass is shifted by 5 cm along x-axis of the Body Frame, i.e., ρ k = 0.05 0 0 m. The initial conditions are presented in Table 2. Here, the satellite free motion is considered, and the initial conditions for translational motion are zero. Table 2. Initial conditions.

Parameter Value
Satellite angular velocity, ω 0.14 0.      As we can see from , the variation of the presented values is lower and more irregular for smaller integration step. So, the numerical model shows the same conservation laws as the mathematical one. Note that in order to correctly simulate the satellite dynamics, it is necessary to consider rather small time step, even for the relatively low RWs angular velocities. In real space applications the RW angular velocity reaches hundreds of radians per second. Additionally, the model includes a number of parameters: inertia tensor, RWs' placement vectors, imbalance and inertia parameters, RWs' angular velocities, as well as control coefficients, which gives at least 24 independent parameters for a satellite with three identical RWs. To have informative enough numerical results, a vast amount of numerical experiments should be carried out. So, the numerical analysis becomes rather computationally expensive. Additionally, the problem of data visualization arises. On the other hand, the analytical estimations on the attitude and stabilization accuracy can be obtained. They allow avoiding extensive numerical simulations during the satellite preliminary design stage. The following sections provide a rather simple end-form expression for the inertial stabilization accuracy.

Effect of Dynamic Imbalances on Inertial Stabilization
In this section, we analyze the effect of the dynamic imbalances in the case of inertial stabilization when the satellite is in the specific attitude with zero angular velocity. Without loss of generality, let the desired attitude be the identity quaternion Q d = 1 0 0 0 T .
All the external torques, as well as RWs' friction, are neglected during the analysis. In addition, we do not consider static imbalance, i.e., RWs' centers of mass are located at their rotation axes.

First-Order Equations of Motion
The dynamic imbalance appears when RW rotation axis is misaligned with the principal axis of inertia, i.e., there are nondiagonal elements in the RW tensor of inertia I k .
The magnitude of these elements d (such that I kj < d for k = j) is usually referred to as dynamic imbalance value. It is rather small in comparison with the RW axial moment of inertia I ax , therefore, we can introduce small parameter ε = d/I ax . The RW tensor of inertia is then represented in the following way: where ε 1, I k is nominal RW tensor of inertia such that RW rotation axis e k is its axis of dynamical symmetry, εδI k is the imbalance additional term. Since I k is axially symmetrical, it does not change in the Body Frame during the RW rotation. Total satellite tensor of inertia in the case of absent static imbalance (i.e., ρ kc = 0) is: Again, here, J is the nominal satellite full tensor of inertia, and ε∑ k δI k corresponds to the small deviations caused by imbalances. Usually, these imbalances are unknown, so simplified equations of motions such as Equation (6) with nominal values of RW and satellite tensors of inertia are used to construct attitude controller.
To analyze the effect of the dynamical imbalances consider equations of motion, Equation (5), taking into account the differences between the nominal and real tensors of inertia. The orbital dynamics are decoupled, in addition, the effect of external torques can be neglected since the time intervals are small. Thus, the system to be considered is: To eliminate a small parameter in the left part of the equations, let us rewrite the first two equations as follows: ). Finally, this system becomes: I 1 e 1 · · · I n e n e T 1 I 1 e T 1 I 1 e 1 · · · 0 . . . The main idea is to represent this system as: where B does not contain a small parameter ε. Detailed derivation is provided in the Appendix B, and here we just present the result: where These equations are used later to analyze the satellite motion in the vicinity of the desired position.

Controller Design
Lyapunov-based attitude controller that ensures the convergence to the necessary [31,32] attitude is: Here, the simplicity of the desired motion (i.e., desired angular velocity and acceleration are equal to zero) is already taken into account. In addition, we consider the case when there are no external disturbances affecting the motion. Note that there is a nominal tensor of inertia in the equation because the attitude control system does not have information about small deviations in the satellite inertia. This is the ideal control torque to be produced by the system of RWs. The necessary RWs' angular accelerations are then defined by: Note that: Ω k e k , I k Ω k = I k Ω k e k Introducing: G = I 1 e 1 . . . I n e n , Ω = Ω 1 . . . Ω n T Equation (9) is rewritten as follows: This system of linear equations allows us to calculate the RWs' accelerations using information about current RWs' angular velocities and the necessary control torque. To ensure controllability of the system, it is necessary to install at least three RWs with noncoplanar rotation axes. In this case, Expression (10) has a unique solution: .
The solution is not unique if there are more than three RWs. In this case, it is reasonable to set an optimization problem: Note that if there are only three RWs, then G T GG T −1 = G −1 , and the solution is the same as in Equation (11). In a simple model RWs' accelerations are: Hence, internal control torques are defined by: If all RWs in the system are identical (at least their nominal axial moments of inertia I k = I), then: G = I e 1 . . . e n = IA and (12) becomes:

Equations of Motion Analysis
Considered controller ensures asymptotic stability of simple system Equation (7), hence we can expect that if ε is sufficiently small, the satellite motion would be in the vicinity of the required one. Therefore, it is possible to linearize the equations of motion. Then, the attitude quaternion is: where ϕ = ϕ 1 ϕ 2 ϕ 3 correspond to three Euler rotation angles (sequence 1-2-3 at angles ϕ 1 , ϕ 2 , ϕ 3 ). From Kinematics (3), we also obtain that the satellite angular velocity in linear approximation is ω = .
ϕ. The solution of System (8) and linearized kinematics is represented in the power series: ω = ω 0 + εω 1 + . . . Ω k = Ω 0 k + εΩ 1 k + . . . = e k Ω 0 k + εΩ 1 k + . . . q ≈ 1 2 ϕ 0 + εϕ 1 + . . . . Therefore, RWs' controls become: Their substitution in the equations of motion and comparison of terms with the same powers of ε gives the system for the undisturbed motion (zero approximation): Zero approximation has an asymptotically stable solution ω 0 = 0, ϕ 0 = 0. This leads to: In addition, from the second equation of zero approximation, we can see that Ω 0 k → const. Therefore, the equations for the first approximation are simplified: Note that δI k is not constant in the Body Frame since the imbalanced RWs rotate. It can be represented in the Body Frame as follows: where: which can be solved analytically. Note that α 0 k is a constant initial phase. The value of the phase depends on the actual transient motion before the satellites settle near the required state. The solution of the homogeneous equation converges to zero asymptotically while the partial solution is: where: The angular velocity for the first approximation is: Hence the attitude stabilization error can be estimated by: For each i-th component, then: This simple estimate might be utilized at the preliminary stage of a spacecraft design to determine the dynamic imbalance requirements for RWs. Each term in Equation (18) is equivalent to: when Ω 0 k is large enough, here, d is the dynamic imbalance magnitude.

Illustrative Example
Assume that each RW tensor of inertia in its own Frame is: where I eq , I ax are the RWs' equatorial and axial moments of inertia (all RW nominal values are supposed to be identical). Note that ∆I RW k is symmetric, and all its components are considered small with respect to the axial moment of inertia.
In RW Frame: Let for all RWs ∆I k 13 = 0, and ∆I k 13 = −d is the dynamic imbalance. Small parameter is introduced as follows: ε = d/I ax Let the nominal RWs' axes of rotation coincide with J principal axes of inertia. Then, in the Body Frame: Matrices that describe the rotation from the RW Frame to the Body Frame are: Here α k = α 0 k + Ω 0 k t, as in the previous section, corresponds to the current RW rotation angle. Finally, the right part of the Equation (15) becomes: In this case, vectors f k and g k in Equation (15) are:

I ax
For the illustration purposes, consider the evolution of the angular velocity after the transient motion. So, the satellite angular velocity is zero, quaternion is identical and almost all initial angular momentum is stored in the RWs. System parameters are presented in Table 1. The control parameters are k ω = 0.01 N·m·s, k a = 0.001 N·m. The results are presented in Figures 6-8. The red curves in figures are the numerical solutions, the blue curves are the approximate solutions in Equation (16), the black horizontal lines are the estimations in Equation (18). In order to test a more realistic scenario, we also include gravity gradient torque into the simulation. This torque is rather small, so at sufficiently small time spans, it would not affect the results.  (16), the black horizontal lines are the estimations in Equation (18). In order to test a more realistic scenario, we also include gravity gradient torque into the simulation. This torque is rather small, so at sufficiently small time spans, it would not affect the results.      First of all, one can see from figures that estimations in Equation (18) are in good accordance with the numerical simulation: numerical results are within the estimation borders, and relative difference between the linearized model of motion and the full one is around 3%. These estimations are of the utmost practical interest since these values show the stabilization accuracy of the satellite. The evolution of the angular velocity components is also close to the numerical solution. However, it should be noted that for the first and second components (Figures 6 and 7), the phase difference increases by the end of the time interval. This is due to the non-uniform evolution of the RWs' angles of rotation. Figure 9 illustrates this effect. First of all, one can see from figures that estimations in Equation (18) are in good accordance with the numerical simulation: numerical results are within the estimation borders, and relative difference between the linearized model of motion and the full one is around 3%. These estimations are of the utmost practical interest since these values show the stabilization accuracy of the satellite. The evolution of the angular velocity components is also close to the numerical solution. However, it should be noted that for the first and second components (Figures 6 and 7), the phase difference increases by the end of the time interval. This is due to the non-uniform evolution of the RWs' angles of rotation. Figure 9 illustrates this effect. is around 3%. These estimations are of the utmost practical interest since these values show the stabilization accuracy of the satellite. The evolution of the angular velocity components is also close to the numerical solution. However, it should be noted that for the first and second components (Figures 6 and 7), the phase difference increases by the end of the time interval. This is due to the non-uniform evolution of the RWs' angles of rotation. Figure 9 illustrates this effect. Figure 9. Phase difference between uniform and numerical RW angle evolution. Figure 9. Phase difference between uniform and numerical RW angle evolution.
In Figure 9, one can see that the phase difference for the third RW is large, which results in the noticeable difference between the numerical and approximate solutions for the first and second angular velocity components.
The following numerical example is considered for the Relation (19) illustration. The first RW angular velocity is taken as Ω 1 = 10 rad/s for the first case, Ω 1 = 50 rad/s for the second one and Ω 2 = Ω 3 = 0 for both cases. Other parameters are the same. Results are presented in Figure 10. In Figure 9, one can see that the phase difference for the third RW is large, which results in the noticeable difference between the numerical and approximate solutions for the first and second angular velocity components.
The following numerical example is considered for the Relation (19) illustration. The first RW angular velocity is taken as 1 10 rad/s Ω = for the first case, 1 50 rad/s Ω = for the second one and 2 3 0 Ω = Ω = for both cases. Other parameters are the same. Results are presented in Figure 10. From Figure 10, one can see that if 1 Ω increases fivefold then the amplitude of angular velocity oscillations also increases fivefold (the angular velocity amplitude increases by the factor of 5.04). So, the Relation (19) is in good accordance with numerical simulation.

Effect of Static Imbalance on Inertial Stabilization
We apply the similar technique to the study of the effect of static imbalance. Its value is usually higher than the one for the dynamic imbalance. In addition, it turns out that the stabilization accuracy in this case depends not only on the pure imbalance parameters.
Using the same approach as in Section 3, the following first order approximation equations are obtained (here k I are supposed to be axially symmetrical): From Figure 10, one can see that if Ω 1 increases fivefold then the amplitude of angular velocity oscillations also increases fivefold (the angular velocity amplitude increases by the factor of 5.04). So, the Relation (19) is in good accordance with numerical simulation.

Effect of Static Imbalance on Inertial Stabilization
We apply the similar technique to the study of the effect of static imbalance. Its value is usually higher than the one for the dynamic imbalance. In addition, it turns out that the stabilization accuracy in this case depends not only on the pure imbalance parameters.
Using the same approach as in Section 3, the following first order approximation equations are obtained (here I k are supposed to be axially symmetrical): Here δρ kc represents the vector from point O k to the RW center of mass and depends on RW rotation angle: where D k is determined by Equation (14). The resulting equations are similar to Equation (15), so the resulting equations of motion are also similar.
Note that attitude stabilization accuracy depends on the RW position with respect to the system center of mass ρ k , i.e., the farther RWs are placed from the system center of mass, the worse stabilization accuracy is. This result is especially useful as it allows us to reduce the effect of vibrations caused by static imbalance at the early stage of satellite design.
Consider  Table 1) for all RWs. All other parameters are the same (without dynamic imbalance). This leads to: In this case vector f k and g k in Equation (15) are: The results of illustrative numerical simulation are presented in Figures 11-13. Black lines are the analytical estimations, red and blue lines are the numerical simulation results.
The example shows that the Expressions (21) are in a good accordance with numerical simulations. As one can see from the close-ups, the difference between analytical estimations and numerical simulation results is rather small and is constrained to 5% in the worst case. The difference can be explained by the gravity gradient torque which is included in the numerical model but is not taken into account in the analytical study.
As one can see, both dynamic and static imbalances lead to the additional terms in the right parts of Equations (15) and (20), which do not depend on satellite state vector components for the first order approximation, so the total stabilization accuracy is the sum of Equations (17) and (21). The results of illustrative numerical simulation are presented in Figures 11-13. Black lines are the analytical estimations, red and blue lines are the numerical simulation results.      The example shows that the Expressions (21) are in a good accordance with numerical simulations. As one can see from the close-ups, the difference between analytical estimations and numerical simulation results is rather small and is constrained to 5% in the worst case. The difference can be explained by the gravity gradient torque which is included in the numerical model but is not taken into account in the analytical study.
As one can see, both dynamic and static imbalances lead to the additional terms in the right parts of Equations (15) and (20), which do not depend on satellite state vector components for the first order approximation, so the total stabilization accuracy is the sum of Equations (17) and (21).

Conclusions
In this paper, the satellite motion analysis is carried out. The model includes RWs' static and dynamic imbalances along with coupling orbital and angular satellite motion with RW rotation. Software implementation of the model is validated using the momentum, angular momentum, and kinetic energy conservation laws. The simulations show that the integration step should be rather small due to the high values of typical RW angular velocities. This fact makes purely numerical analysis difficult. In order to solve this problem, the analytical approximations for the satellite stabilization accuracy are obtained in closed form for the static and dynamic imbalances' presence in the inertial stabilization case. The comparison of the numerical simulation and approximate solution shows that they are in a good accordance (relative error is about several percent). The explicit expressions can easily be implemented and are useful during the preliminary satellite design stage.
The estimations of the attitude accuracy are obtained for the case of the satellite's inertial stabilization. The case of orbital stabilization is the main goal of future research.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In order to derive equations of motion, the general equation of dynamics is used. Terms for each independent virtual displacements are: We start with the equations of motion for point O (first equation). Since: ..
it can be rewritten as follows: where: Ω k + ω × Ω k × ρ kc + ∑ k m k (ω + Ω k ) × (ω + Ω k ) × ρ kc = F s + ∑ k F k . Introducing: and utilizing properties: equations of hull angular motion are derived from the second Equation of (A1): Under Constraint (A2) and using total satellite tensor of inertia: From the third equation of (A1), we derive: e T k m k ρ kc × ..
After simplification: e T k m k ρ kc × ..

Appendix B
We consider the system of equations: I k e k = I k e k , e T k I k e k = I k , e T k δI k e k = δI k , and want to eliminate small parameter ε near the highest order derivatives . ω, . Ω k . Since Ω k = e k Ω k , from the second Equation: . Ω k = b k +εδb k −e T k ( I k +εδI k ) . ω ( I k +εδI k ) = 1