State-Space Motor Modeling

The questions below are due on Thursday April 25, 2024; 10:30:00 PM.

You are not logged in.
Please Log In for full access to the web site.
Note that this link will take you to an external site (https://shimmer.mit.edu) to authenticate, and then you will be redirected back to this page.

1) An L-Motor-Driven Arm

In this problem, we will consider an example that is similar to the propeller-levitated arm, a lego-motor-driven arm (for brievity, an Lmotor-driven arm). The similarities will make it easier for us to understand the structure of this new control problem, but differences will allow us to raise additional issues, while also giving us a chance to practice with the state-space concepts.

In earlier versions of 6.302, we rotated a heavier arm using a geared Lego motor, shown in the figure below. The figure also shows that we used two rubber bands to create an angle dependent torque on the motor (both to assist the motor when holding heavy loads, and to eliminate gear "backlash"). That angle dependent force will have a number of interesting consequences, ones that we hope will help you develop a better understanding of controller design.

The L-motor is a geared brushed motor, just like the propeller motor. By changing the motor voltage, v, we change the motor current, i_m, which is proportional to the torque generated by the motor. The net arm torque, a sum of motor-, gravity-, and elastic-generated torques, results in an arm rotation characterized angular acceleration, \alpha, and angular velocity, \omega, and rotation angle \theta, as shown in the figure below.

We can sense the Lmotor-driven arm angle using the same kind of angle sensor we use for the propeller arm, by linking the shaft of the angle sensor to the shaft of the motor (using heat-shrink tubing actually). As an example, using our angle sensor, we could easily meaure the \theta in the figure below as approximately \frac{-\pi}{6}.

In fact, we can reuse most of our propeller-levitated arm hardware in an Lmotor-driven arm controller, since they both use brushed DC motors. We can use the same angle sensor, and we can drive the Lmotor using the same H-bridge (the Pololu board) with almost the same PWM voltages (though lego motors are designed to work with lower-frequency PWM signals). In fact, we could even use exactly the same PID controller, though we might find that good K_i, K_d, and K_p values for the Lmotor-driven arm will be different from good values for the propeller arm.

The issues in designing state-space controllers are a little different, however, as we shall see.

2) Generating a State-Space model

We will need a few physical relationships in order to derive a state-space model, but thankfully, the Lmotor-driven arm model does not require linearization.

2.1) The physical relationships

As in the propeller-levitated arm, the Lmotor-driven arm angle, \theta(t), is related to its angular velocity, \omega(t), through differentiation,

\frac{d \theta}{dt} = \omega(t),

and of course, the angular velocity of the arm is related to arm angular acceleration through differentiation,

\frac{d \omega(t)}{dt} = \alpha(t).

The product of the arm's moment of inertia, J, and its angular acceleration, \alpha, is equal to the torque on the arm \tau,

J \alpha(t)=\tau(t) = \tau_m(t) + \tau_b(t) + \tau_f(t),

where \tau_m is the motor-generated torque, \tau_b is the net torque generated by the rubber bands, and \tau_f is the torque due to friction. We are going to assume the arm is rotating parallel to the ground, so that we can ignore gravity as a source of torque on the arm.

As we have seen before, motor torque, \tau_m, is proportional to motor current,

\tau_m = K_m i_m,

where i_m is the motor current and K_m is the current to torque conversion factor. We can model the net torque generated by the two rubber bands as proportional to arm angle, as long as the bands have been placed so that their forces balance at \theta =0. If so, then

\tau_b = K_b \theta,

where \tau_b is the net elastic torque and K_b is the proportionality constant. Finally, the friction torque can be modeled as proportional to angular velocity,

\tau_f = -K_f \omega,

where K_f is the friction coefficient, and the negative sign is because the friction force always acts to oppose motion.

Since the Lmotor is a brushed motor, we can model it in the same way we modeled the propeller motor. That is, we can use a circuit model to relate the voltage across the motor to the current through the motor, and after scaling by K_m, relate the motor voltage to the motor torque. If we model the motor as a series combination of a motor coil resistor, R_m, a motor coil inductor, L_m, and a rotational-velocity dependent voltage source, v_{emf}, the current through the motor, i_m, can be related to the drive voltage, v, by

v(t) = R_m i_m(t) + L_m\frac{di_m(t)}{dt} + v_{emf}(t).

The motor's inductance and resistance are much larger than in the propeller motor, so the above relation can not be simplified. The v_{emf} is proportional to the motor speed, \omega_m, a fact we made use of, but did not investigate, in the SPINUP sketch used to calibrate our propeller motor.

v_{emf}(t) = K_e\omega(t).

2.2) Values for the physical parameters.

(http://nxt-unroller.blogspot.com/2015/03/mathematical-model-of-lego-ev3-motor.html)

L_m = 0.005 is the motor inductance in Henries.

J = 0.0015 is the arm-plus-motor moment of inertia.

R_m = 7.0 is the motor's series resistance in ohms.

K_e = 0.46 is the back-emf per radian/sec of motor rotational velocity.

K_m = 0.3 is the torque per ampere of current.

K_f = 0.00073 is the torque due to friction per radian/second.

K_b = 0.0 \;or\; -0.1 is the elastic torque per radian (due to the rubber bands) and is negative because the elastic was used to oppose arm motion.

2.3) Generating the State-Space Matrices

We can describe our Lmotor-driven arm, our "plant", as a single-input single -output (SISO) state-space system with the motor voltage as the input and the arm angle as the output. The generic form of the state-space model using \textbf{E}, \textbf{A}, \textbf{B}, and \textbf{C} matrices is

\textbf{E} \frac{d}{dt} \textbf{x} = \textbf{A}\textbf{x}(t) + \textbf{B} u(t)

and

y(t) = \textbf{C} \textbf{x}(t),

where the states in the system (there are three) are the elements of the three-length vector, \textbf{x}, u is the motor voltage, and y is arm angle.

There are an infinite number of state-space models that can be generated from our physical model of the Lmotor-driven arm, even though we have specified the input and output. To make it easier to grade, we will constrain which element of the vector \textbf{x} goes with which physical state by fixing the input and output matrices, \textbf{B} and \textbf{C},

\textbf{B} = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}

\textbf{C} = \begin{bmatrix} 0&0 & 1\end{bmatrix}.

Unfortunately, there are so many zeros in \textbf{B} and \textbf{C} that the state-space system is still not unique, so we will also constrain the \textbf{E} to be diagonal, with diagonal entries L_m, J, and 1.

What is the state-space matrix \textbf{A} for our Lmotor driven arm, assuming that there are rubber bands? Please enter the \textbf{A} matrix below, and remember to use Python form, so, for example, the \textbf{E} matrix would be specified as

`[[5.0e-3, 0, 0], [0, 1.5e-3,0],[0, 0, 1]]'

3) Poles and Eigenvalues.

In the SISO state-space system, the transfer function from plant input u to plant output y is given by

H_{plant}(s) = \textbf{C} \left(s\textbf{E} - \textbf{A}\right)^{-1} \textbf{B},

but how does that formula relate to the transfer functions from classical control. For example, how do we relate the above matrix formula to an L-zero, N-pole transfer function,

H(s) = \frac{(s - z_1)(s-z_2)...(s-z_{L-1})(s-z_L)}{(s - p_1)(s-p_2)...(s-p_{N-1})(s-p_N)},

which, if the poles are distinct, can be expanded in to

H(s) = \frac{\alpha_1}{s - p_1} + \frac{\alpha_2}{s - p_2} + ... + \frac{\alpha_N}{s - p_N}.

Let us start with the \textbf{E} = \textbf{I} (the identity matrix) case. In this case, the eigenvalues of \textbf{A} are the poles of H(s) Recall that if \lambda_i is an eigenvalue of \textbf{A}, and \textbf{V}_i is the associated eigenvector, then

\textbf{A} \textbf{V}_i = \lambda_i \textbf{V}_i.

If there are no repeated eigenvalues, that is the eigenvalues are distinct, then

\textbf{A} \begin{bmatrix} \textbf{V}_1 \textbf{V}_2 ... \textbf{V}_N \end{bmatrix} = \begin{bmatrix} \textbf{V}_1 \textbf{V}_2 ... \textbf{V}_N \end{bmatrix} \begin{bmatrix} \lambda_1 & 0 & ... & ... & 0\\ 0 & \lambda_2 & 0 & ... & 0 \\ ... & ... & ... & ... & ...\\ 0 & ... & ... & 0 & \lambda_N \end{bmatrix},

\textbf{A} \textbf{V} = \textbf{V} \Lambda

where \textbf{V} is a matrix of eigenvectors. The columns of \textbf{V} are the eigenvectors of \textbf{A}, and are guaranteed to be linearly independent if \textbf{A}'s eigenvalues are distinct, and therefore \textbf{V} will be invertible. Then using the inverse of \textbf{V},

\textbf{V}^{-1} \textbf{A} \textbf{V} = \Lambda,

and we say that the eigenvectors diagonalize \textbf{A}.

Consider a change of variables, in which the vector \textbf{x} is written as a weighted combination of the eigenvectors of \textbf{A} (also known as \textbf{x}'s spectral decomposition),

\textbf{x}(t) = \textbf{V} \textbf{w}(t)

where \textbf{w} is the vector of eigenvector weights. Using this change of variables in the state-space model,

\frac{d}{dt} \textbf{V} \textbf{w}(t) = \textbf{A}\textbf{V} \textbf{w}(t) + \textbf{B} u(t),

and

y(t) = \textbf{C} \textbf{V} \textbf{w}(t).

If we then multiply by \textbf{V}^{-1}, and define \tilde{\textbf{B}} and \tilde{\textbf{C}} in the new variables,

\frac{d}{dt} \textbf{w}(t) = \textbf{V}^{-1} \textbf{A}\textbf{V} \textbf{w}(t) + \textbf{V}^{-1}\textbf{B} u(t) = \begin{bmatrix} \lambda_1 & 0 & ... & ... & 0\\ 0 & \lambda_2 & 0 & ... & 0 \\ ... & ... & ... & ... & ...\\ 0 & ... & ... & 0 & \lambda_N \end{bmatrix} \textbf{w}(t) + \tilde{\textbf{B}} u(t) = \Lambda \textbf{w}(t) + \tilde{\textbf{B}} u(t)

and

y(t) = \textbf{C} \textbf{V} \textbf{w}(t) = \tilde{\textbf{C}} \textbf{w}(t).

where we use the tilde to denote matrices after the change of variables.

Rewriting

H_{plant}(s) = \textbf{C} \left(s\textbf{I} - \textbf{A}\right)^{-1} \textbf{B} = \tilde{\textbf{C}} \begin{bmatrix} \frac{1}{s - \lambda_1} & 0 & ... & ... & 0\\ 0 & \frac{1}{s - \lambda_2} & 0 & ... & 0 \\ ... & ... & ... & ... & ...\\ 0 & ... & ... & 0 & \frac{1}{s - \lambda_N} \end{bmatrix} \tilde{\textbf{B}}.

If we write out the matrix products and match it to the partial fraction expansion from above,

H(s) = \frac{\tilde{\textbf{C}}_1 \tilde{\textbf{B}}_1}{s - \lambda_1} + \frac{\tilde{\textbf{C}}_2 \tilde{\textbf{B}}_2}{s - \lambda_2} + ... + \frac{\tilde{\textbf{C}}_N \tilde{\textbf{B}}_N}{s - \lambda_N} = \frac{\alpha_1}{s - p_1} + \frac{\alpha_2}{s - p_2} + ... + \frac{\alpha_N}{s - p_N},

we can see that the pole p_i equals the eigenvalue \lambda_i, and the coefficient in the partial fraction expansion, \alpha_i equals the product \tilde{\textbf{C}}_i \tilde{\textbf{B}}_i.

Now consider the case where \textbf{E} is not equal to \textbf{I}, in which case we will make use of generalized eigenvectors and eigenvalues. We say that \lambda_i is a generalized eigenvalue and \textbf{V}_i is the associated generalized eigenvector of the matrix pair \textbf{A},\textbf{E} if

\textbf{A} \textbf{V}_i = \lambda_i \textbf{E} \textbf{V}_i.

and if there are N eigenvalues and eigenvectors, then

\textbf{A} \textbf{V} = \textbf{E} \textbf{V} \Lambda

where \Lambda is the diagonal matrix of eigenvalues and \textbf{V} is the matrix of eigenvectors. The generalized eigenvector and eigenvalue matrices can be computed with matlab's 'eig' command using two input arguments, as in

[V,Lambda] = eig(A,E)

If \textbf{E} matrix is invertible, and the generalized eigenvalues are distinct, we can still use a change of variables to convert the state-space model to

\frac{d}{dt} \textbf{w}(t) = \begin{bmatrix} \lambda_1 & 0 & ... & ... & 0\\ 0 & \lambda_2 & 0 & ... & 0 \\ ... & ... & ... & ... & ...\\ 0 & ... & ... & 0 & \lambda_N \end{bmatrix} \textbf{w}(t) + \tilde{\textbf{B}} u(t) = \Lambda \textbf{w}(t) + \tilde{\textbf{B}} u(t)

and

y(t) = \tilde{\textbf{C}} \textbf{w}(t).

3.1) Poles and transfer functions Beyond E = Identity

Assuming \textbf{E} is invertible (but is not necessarily the identity) and that the generalized eigenvalues are distinct, answer the following questions in terms of E, A, B, C, V and L (for \Lambda), and their inverses (e.g. E**-1)

The diagonals of L are the eigenvalues of

In terms of E, V and B, and their inverses, Btilde =

In terms of E, V, and C, and their inverses, Ctilde =

When using state-feedback control

u(t) = K_r r(t) - \textbf{K} \textbf{x}(t)

where r is the reference input. The closed-loop transfer function from reference input to plant output is given by

H_{closedLoop}(s) = \textbf{C} \left(s\textbf{E} - \left(\textbf{A}-\textbf{B}\textbf{K}\right) \right)^{-1} \textbf{B} K_r = \textbf{C} \left(s\textbf{I} - \textbf{E}^{-1}\left(\textbf{A}-\textbf{B}\textbf{K}\right) \right)^{-1} \textbf{E}^{-1}\textbf{B} K_r.

The poles of the closed-loop transfer function are equal to the eigenvalues of what matrix (use A**-1 to denote A^{-1})?

Poles of the closed-loop transfer function are the eigenvalues of what matrix?

3.2) Poles and transfer functions for the Lmotor-driven arm

Now use the above formulas (and matlab) to compute the poles and partial fraction expansion representation of the transfer function for the open-loop Lmotor-driven arm plant,

H(s) = \frac{\alpha_1}{s - p_1} + \frac{\alpha_2}{s - p_2} + \frac{\alpha_3}{s - p_3},

which can be reorganized as

H(s) = \frac{\beta_1}{\frac{s}{p_1} - 1} + \frac{\beta_2}{\frac{s}{p_2} - 1} +\frac{\beta_3}{\frac{s}{p_3} - 1},

where \beta_i = \frac{\alpha_i}{p_i}.

If Kb=0 (no rubber bands), what are the poles of the Lmotor-driven-arm plant (enter in python list format, e.g. [4.0,5.0+3.0j,5.0-3.0j])?

If Kb=0 (no rubber bands), one of the poles in the Lmotor-driven arm plant is zero, and therefore its associated beta is infinite. What are the two remaining finite beta's (enter in python list format, e.g. [4.0,5.0+3.0j,5.0-3.0j])?

If Kb= -0.1 (with rubber bands), what are the poles of the Lmotor-driven arm plant (enter in python list format, e.g. [4.0,5.0+3.0j,5.0-3.0j])?

If Kb = -0.1 (with rubber bands), what are the beta's defined above for the Lmotor-driven-arm plant (enter in python list format, e.g. [4.0,5.0+3.0j,5.0-3.0j])?

Note that one pole is much larger than the other two, and its associated \beta much smaller. That suggests H(s) could be well approximated without that pole, an idea we return to below.

4) The elastic torque and K_r

Since you completed the above problem, you know that the states for the system are i_m, \omega, and \theta, in that order. If you want to design a state-feedback controller, then you have to determine K_r and \textbf{K} for

u = K_r y_d - \textbf{K} \textbf{x}

or in this case

v = K_r \theta_d - \left( \textbf{K}_{1} i_m + \textbf{K}_{2} \omega + \textbf{K}_3 \theta \right),

where \theta_d is the desired arm angle.

Recall that given a vector state feedback gains, \textbf{K}, we derived a formula for K_r based on matching the steady-state for the reference input (if it has one) to the steady-state of the output. But for the propeller arm examples we have done in lab, the formula for K_r,

K_r = \frac{-1}{\left(\textbf{C} (\textbf{A} - \textbf{B}\textbf{K})^{-1} \textbf{B}\right)},

always resulted in a K_r that was always equal to one of the state-feedback gains. What we saw in lab is not always the case, but often is, and with a little linear algebra trick we will see why.

When we create SISO state-space models for systems of interest, we often find it convenient to pick representations in which only one of the elements of \textbf{B} and one of the elements of \textbf{C} are nonzero, and those nonzero elements are equal to one. We did exactly that for the Lmotor-driven arm model, by setting

\textbf{B} = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}, \;\;\; \textbf{C} = \begin{bmatrix} 0&0 & 1\end{bmatrix}.

For our specific \textbf{B} and \textbf{C}, the value of K_r depends on only one element of the inverse of (\textbf{A} - \textbf{B}\textbf{K})^{-1}. For what i,j is

K_r = \frac{-1}{(\textbf{A} - \textbf{B}\textbf{K})^{-1}_{i,j}}

Enter your answer in python format. So if your answer is i=2 and j=3, enter

[2,3]

[i,j]=

When calculating the inverse of a matrix \textbf{M}, it is sometimes helpful to recall that if the inverse of the matrix exists, then

\textbf{M} \textbf{M}^{-1} = \begin{bmatrix} 1 & 0 & ... & ... & 0\\ 0 & 1 & 0 & ... & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & ... & 0 & 1 & 0\\ 0 & ... & 0 & 0 & 1 \end{bmatrix} = \textbf{I}

In particular, if we are interested in computing the first column of M^{-1}, all we have to do is find a vector \textbf{w} for which

\textbf{M} \textbf{w} = \begin{bmatrix} 1 \\ 0 \\ \vdots\\ 0 \end{bmatrix}.

This column at a time approach to finding inverses is particularly helpful when the nonzero entries of \textbf{M} have a certain structure. For example, suppose \textbf{M} has an inverse, and in addition, the only nonzero element of the 5^{th} column of \textbf{M} is M_{1,5} = 25. For this example, the first column of \textbf{M}^{-1} will have only one nonzero (do you see why?). In what row is that nonzero, and what is its value?

row =

value =

The above observation will help you understand when K_r will be equal to one of the state gains. In particular, suppose you remove the elastic bands from the motor. Then K_b = 0 in the model, and in the open loop case, nothing depends directly on arm angle. For this K_b = 0 case, give an expression for K_r in terms of the three state gains (in answering the question, denote the three state gains as 'K1', 'K2' and 'K3').

Kr=

Suppose the elastic bands are back on the motor, so that K_b = -0.1, and the arm torque is still dependent on state gains are

\textbf{K} = \begin{bmatrix} 1&2 & 10\end{bmatrix},

what is the numerical value of K_r?

Kr=

The model for the propeller-levitated arm we used in lab is like the Lmotor-driven arm without elastic bands. Nothing in the propeller-arm model depended directly on arm angle, and as a result, K_r was always equal to K_{angle}.