To solve this problem using Matlab, copy the following commands into an new m-file: m=1000; b=50; u=500; num=[1]; den=[m b]; These commands will later be used to find the open-loop response of the system to a step input. But before getting into that, let's take a look at another representation, the state-space. 2. State-Space We can rewrite the first-order modeling equation (1) as the state-space model. To use Matlab to solve this problem, create an new m-file and copy the following commands: m b u A B C D = = = = = = = 1000; 50; 500; [-b/m]; [1/m]; [1]; 0; Note: It is possible to convert from the state-space representation to the transfer function or vise versa using Matlab. To learn more about the conversion, click Conversion Open-loop response Now let's see how the open-loop system responds to a step input. Add the following command onto the end of the m-file written for the tranfer function (the m-file with num and den matrices) and run it in the Matlab command window: step (u*num,den) You should get the following plot: To use the m-file written for the state-space (the m-file with A, B, C, D matrices), add the following command at the end of the m-file and run it in the Matlab command window: step (A,u*B,C,D) You should get the same plot as the one shown above. From the plot, we see that the vehicle takes more than 100 seconds to reach the steady-state speed of 10 m/s. This does not satisfy our rise time criterion of less than 5 seconds.
law Closed-loop transfer function To solve this problem, a unity feedback controller will be added to improve the system performance. The figure shown below is the block diagram of a typical unity feedback system. The transfer function in the plant is the transfer function derived above {Y(s)/U(s)=1/ms+b}. The controller will to be designed to satisfy all design criteria. Four different methods to design the controller are listed at the bottom of this page. You may choose on PID, Root-locus, Frequency response, or State-space. Example: DC Motor Speed Modeling Physical setup and system equations Photo courtesy: Pope Electric Motors Pty Limited A common actuator in control systems is the DC motor. It directly provides rotary motion and, coupled with wheels or drums and cables, can provide transitional motion. The electric circuit of the armature and the free body diagram of the rotor are shown in the following figure: For this example, we will assume the following values for the physical parameters. These values were derived by experiment from an actual motor in Carnegie Mellon's undergraduate controls lab. * moment of inertia of the rotor (J) = 0.01 kg.m^2/s^2 * damping ratio of the mechanical system (b) = 0.1 Nms * electromotive force constant (K=Ke=Kt) = 0.01 Nm/Amp * electric resistance (R) = 1 ohm * electric inductance (L) = 0.5 H * input (V): Source Voltage * output (theta): position of shaft * The rotor and shaft are assumed to be rigid The motor torque, T, is related to the armature current, i, by a constant factor Kt. The back emf, e, is related to the rotational velocity by the following equations: In SI units (which we will use), Kt (armature constant) is equal to Ke (motor constant). From the figure above we can write the following equations based on Newton's combined with Kirchhoff's law: 1. Transfer Function Using Laplace Transforms, the above modeling equations can be expressed in terms of s.
• • • By eliminating I(s) we can get the following open-loop transfer function, where the rotational speed is the output and the voltage is the input. 2. State-Space In the state-space form, the equations above can be expressed by choosing the rotational speed and electric current as the state variables and the voltage as an input. The output is chosen to be the rotational speed. Design requirements First, our uncompensated motor can only rotate at 0.1 rad/sec with an input voltage of 1 Volt (this will be demonstrated later when the open-loop response is simulated). Since the most basic requirement of a motor is that it should rotate at the desired speed, the steady-state error of the motor speed should be less than 1%. The other performance requirement is that the motor must accelerate to its steady-state speed as soon as it turns on. In this case, we want it to have a settling time of 2 seconds. Since a speed faster than the reference may damage the equipment, we want to have an overshoot of less than 5%. If we simulate the reference input (r) by an unit step input, then the motor speed output should have: Settling time less than 2 seconds Overshoot less than 5% Steady-state error less than 1% Matlab representation and open-loop response 1. Transfer Function We can represent the above transfer function into Matlab by defining the numerator and denominator matrices as follows: Create a new m-file and enter the following commands: J=0.01; b=0.1; K=0.01; R=1; L=0.5; num=K; den=[(J*L) ((J*R)+(L*b)) ((b*R)+K^2)]; Now let's see how the original open-loop system performs. Add the following commands onto the end of the m- file and run it in the Matlab command window: step(num,den,0:0.1:3) title('Step Response for the Open Loop System') You should get the following plot:
From the plot we see that when 1 volt is applied to the system, the motor can only achieve a maximum speed of 0.1 rad/sec, ten times smaller than our desired speed. Also, it takes the motor 3 seconds to reach its steady-state speed; this does not satisfy our 2 seconds settling time criterion. 2. State-Space We can also represent the system using the state-space equations. Try the following commands in a new m-file. J=0.01; b=0.1; K=0.01; R=1; L=0.5; A=[-b/J -K/L B=[0 K/J -R/L]; 1/L]; C=[1 0]; D=0; step(A, B, C, D) Run this m-file in the Matlab command window, and you should get the same output as the one shown above. Example: Modeling DC Motor Position Physical Setup A common actuator in control systems is the DC motor. It directly provides rotary motion and, coupled with wheels or drums and cables, can provide transitional motion. The electric circuit of the armature and the free body diagram of the rotor are shown in the following figure: For this example, we will assume the following values for the physical parameters. These values were derived by experiment from an actual motor in Carnegie Mellon's undergraduate controls lab. * moment of inertia of the rotor (J) = 3.2284E-6 kg.m^2/s^2 * damping ratio of the mechanical system (b) = 3.5077E-6 Nms * electromotive force constant (K=Ke=Kt) = 0.0274 Nm/Amp * electric resistance (R) = 4 ohm * electric inductance (L) = 2.75E-6 H * input (V): Source Voltage * output (theta): position of shaft * The rotor and shaft are assumed to be rigid System Equations The motor torque, T, is related to the armature current, i, by a constant factor Kt. The back emf, e, is related to the rotational velocity by the following equations:
• • • • In SI units (which we will use), Kt (armature constant) is equal to Ke (motor constant). From the figure above we can write the following equations based on Newton's law combined with Kirchhoff's law: 1. Transfer Function Using Laplace Transforms the above equations can be expressed in terms of s. By eliminating I(s) we can get the following transfer function, where the rotating speed is the output and the voltage is an input. However during this example we will be looking at the position, as being the output. We can obtain the position by integrating Theta Dot, therefore we just need to divide the transfer function by s. 2. State Space These equations can also be represented in state- space form. If we choose motor position, motor speed, and armature current as our state variables, we can write the equations as follows: Design requirements We will want to be able to position the motor very precisely, thus the steady-state error of the motor position should be zero. We will also want the steady-state error due to a disturbance, to be zero as well. The other performance requirement is that the motor reaches its final position very quickly. In this case, we want it to have a settling time of 40ms. We also want to have an overshoot smaller than 16%. If we simulate the reference input (R) by a unit step input, then the motor speed output should have: Settling time less than 40 milliseconds Overshoot less than 16% No steady-state error No steady-state error due to a disturbance
Matlab representation and open-loop response 1. Transfer Function We can put the transfer function into Matlab by defining the numerator and denominator as vectors: Create a new m-file and enter the following commands: J=3.2284E-6; b=3.5077E-6; K=0.0274; R=4; L=2.75E-6; num=K; den=[(J*L) ((J*R)+(L*b)) ((b*R) +K^2) 0]; Now let's see how the original open-loop system performs. Add the following command onto the end of the m-file and run it in the Matlab command window: step(num,den,0:0.001:0.2) You should get the following plot: From the plot we see that when 1 volt is applied to the system, the motor position changes by 6 radians, six times greater than our desired position. For a 1 volt step input the motor should spin through 1 radian. Also, the motor doesn't reach a steady state which does not satisfy our design criteria 2. State Space We can put the state space equations into Matlab by defining the system's matrices as follows: J=3.2284E-6; b=3.5077E-6; K=0.0274; R=4; L=2.75E-6; A=[0 1 0 0 -b/J K/J 0 -K/L -R/L]; B=[0 ; 0 ; 1/L]; C=[1 0 0]; D=[0]; The step response is obtained using the command step(A,B,C,D) Unfortunately, Matlab responds with Warning: Divide by zero ??? Index exceeds matrix dimensions. Error in ==> /usr/local/lib/matlab/toolbox/control/step.m On line 84 ==> dt = t(2)-t(1); There are numerical scaling problems with this representation of the dynamic equations. To fix the problem, we scale time by tscale = 1000. Now the output time will be in milliseconds rather than in seconds. The equations are given by tscale = 1000; J=3.2284E-6*tscale^2;
M m b l I F x theta • • b=3.5077E-6*tscale; K=0.0274*tscale; R=4*tscale; L=2.75E-6*tscale^2; A=[0 1 0 0 -b/J K/J 0 -K/L -R/L]; B=[0 ; 0 ; 1/L]; C=[1 0 0]; D=[0]; The output appears the same as when obtained through the transfer function, but the time vector must be divided by tscale. [y,x,t]=step(A,B,C,D); plot(t/tscale,y) ylabel('Amplitude') xlabel('Time (sec)') Example: Modeling an Inverted Pendulum Problem setup and design requirements The cart with an inverted pendulum, shown below, is "bumped" with an impulse force, F. Determine the dynamic equations of motion for the system, and linearize about the pendulum's angle, theta = Pi (in other words, assume that pendulum does not move more than a few degrees away from the vertical, chosen to be at an angle of Pi). Find a controller to satisfy all of the design requirements given below. For this example, let's assume that mass of the cart 0.5 kg mass of the pendulum 0.5 kg friction of the cart 0.1 N/m/sec length to pendulum center of mass 0.3 m inertia of the pendulum 0.006 kg*m^2 force applied to the cart cart position coordinate pendulum angle from vertical For the PID, root locus, and frequency response sections of this problem we will be only interested in the control of the pendulums position. This is because the techniques used in these tutorials can only be applied for a single-input-single-output (SISO) system. Therefore, none of the design criteria deal with the cart's position. For these sections we will assume that the system starts at equilibrium, and experiences an impulse force of 1N. The pendulum should return to its upright position within 5 seconds, and never move more than 0.05 radians away from the vertical. The design requirements for this system are: Settling time of less than 5 seconds. Pendulum angle never more than 0.05 radians from the vertical. However, with the state-space method we are more readily able to deal with a multi-output system. Therefore, for this section of the Inverted Pendulum example we will attempt to control both the pendulum's angle and the cart's position. To make the design more challenging we will be applying a step input to the cart. The cart should achieve it's desired position within 5 seconds and have a rise time under 0.5 seconds. We will also limit the pendulum's overshoot to 20 degrees (0.35 radians), and it should also settle in under 5 seconds. The design requirements for the Inverted Pendulum state-space example are:
• • • Settling time for x and theta of less than 5 seconds. Rise time for x of less than 0.5 seconds. Overshoot of theta less than 20 degrees (0.35 radians). Force analysis and system equations Below are the two Free Body Diagrams of the system. Summing the forces in the Free Body Diagram of the cart in the horizontal direction, you get the following equation of motion: Note that you could also sum the forces in the vertical direction, but no useful information would be gained. Summing the forces in the Free Body Diagram of the pendulum in the horizontal direction, you can get an equation for N: If you substitute this equation into the first equation, you get the first equation of motion for this system: (1) To get the second equation of motion, sum the forces perpendicular to the pendulum. Solving the system along this axis ends up saving you a lot of algebra. You should get the following equation: To get rid of the P and N terms in the equation above, sum the moments around the centroid of the pendulum to get the following equation: Combining these last two equations, you get the second dynamic equation: (2) Since Matlab can only work with linear functions, this set of equations should be linearized about theta = Pi. Assume that theta = Pi + ø (ø represents a small angle from the vertical upward direction). Therefore, cos(theta) = -1, sin(theta) = -ø, and (d(theta)/dt)^2 = 0. After linearization the two equations of motion become (where u represents the input): 1. Transfer Function
is: To obtain the transfer function of the linearized system equations analytically, we must first take the Laplace transform of the system equations. The Laplace transforms are: NOTE: When finding the transfer function initial conditions are assumed to be zero. Since we will be looking at the angle Phi as the output of interest, solve the first equation for X(s), then substituting into the second equation: Re-arranging, the transfer function where, From the transfer function above it can be seen that there is both a pole and a zero at the origin. These can be canceled and the transfer function becomes: 2. State-Space After a little algebra, the linearized system equations equations can also be represented in state-space form:
0 The C matrix is 2 by 4, because both the cart's position and the pendulum's position are part of the output. For the state-space design problem we will be controlling a multi-output system so we will be observing the cart's position from the first row of output and the pendulum's with the second row. Matlab representation and the open-loop response 1. Transfer Function The transfer function found from the Laplace transforms can be set up using Matlab by inputting the numerator and denominator as vectors. Create an m-file and copy the following text to model the transfer function: M m b i g l = = = = = = .5; 0.2; 0.1; 0.006; 9.8; 0.3; q = (M+m)*(i+m*l^2)-(m*l)^2; %simplifies input num = [m*l/q 0] den = [1 b*(i+m*l^2)/q Your output should be: num = 4.5455 den = -(M+m)*m*g*l/q -b*m*g*l/q] 1.0000 0.1818 -31.1818 -4.4545 To observe the system's velocity response to an impulse force applied to the cart add the following lines at the end of your m-file: t=0:0.01:5; impulse(num,den,t) axis([0 1 0 60]) Note: Matlab commands from the control system toolbox are highlighted in red. You should get the following velocity response plot: As you can see from the plot, the response is entirely unsatisfactory. It is not stable in open loop. You can change the axis to see more of the response if you need to convince yourself that the system is unstable.
0; 1; 0 0 1. State-Space Below, we show how the problem would be set up using Matlab for the state-space model. If you copy the following text into a m-file (or into a '.m' file located in the same directory as Matlab) and run it, Matlab will give you the A, B, C, and D matrices for the state-space model and a plot of the response of the cart's position and pendulum angle to a step input of 0.2 m applied to the cart. M m b i g l = = = = = = .5; 0.2; 0.1; 0.006; 9.8; 0.3; p = i*(M+m)+M*m*l^2; %denominator for the A and B matricies A = [0 1 0 0; 0 -(i+m*l^2)*b/p (m^2*g*l^2)/p 0 0 0 B = [ 0 -(m*l*b)/p 0; m*g*l*(M+m)/p 0] (i+m*l^2)/p; 0; m*l/p] C = [1 0 0 0; 0 0 1 0] D = [0; 0] T=0:0.05:10; U=0.2*ones(size(T)); [Y,X]=lsim(A,B,C,D,U,T); plot(T,Y) axis([0 2 0 100]) You should see the following output after running the m-file: A = 0 0 1.0000 -0.1818 0 -0.4545 0 2.6727 0 31.1818 0 0 1.0000 0 B = 0 1.8182 0 4.5455 C = 1 0 0 0 0 1 0 0 D = 0 0 The blue line represents the cart's position and the green line represents the pendulum's angle. It is obvious from this plot and the one above that some sort of control will have to be designed to improve the dynamics of the system. Four example controllers are included with these tutorials: PID, root locus, frequency response, and state space. Select from below the one you would like to use.
any Note: The solutions shown in the PID, root locus and frequency response examples may not yield a workable controller for the inverted pendulum problem. As stated previously, when we put this problem into the single- input, single-output framework, we ignored the x position of the cart. The pendulum can be stabilized in an inverted position if the x position is constant or if the cart moves at a constant velocity (no acceleration). Where possible in these examples, we will show what happens to the cart's position when our controller is implemented on the system. We emphasize that the purpose of these examples is to demonstrate design and analysis techniques using Matlab; not to actually control an inverted pendulum. Example: Modeling a Pitch Controller Photo courtesy: Boeing Physical setup and system equations Design requirements Transfer function and state-space Matlab representation and open-loop response Closed-loop transfer function Physical setup and system equations The equations governing the motion of an aircraft are a very complicated set of six non-linear coupled differential equations. However, under certain assumptions, they can be decoupled and linearized into the longitudinal and lateral equations. Pitch control is a longitudinal problem, and in this example, we will design an autopilot that controls the pitch of an aircraft. The basic coordinate axes and forces acting on an aircraft are shown in the figure below: Assume that the aircraft is in steady-cruise at velocity; thus, the thrust and drag cancel out and balance out each other. Also, assume that change does not change the speed of an aircraft under circumstance (unrealistic but simplifies the constant altitude and the lift and weight in pitch angle problem a bit). Under these assumptions, the longitudinal equations of motion of an aircraft can be written as: (1) Please refer to any aircraft-related textbooks for the explanation of how to derive these equations. Also, click Variables to see what each variable represents. For this system, the input will be the elevator deflection angle, and the output will be the pitch angle. Design requirements The next step is to set some design criteria. We want to design a feedback controller so that the output has an overshoot of less than 10%, rise time of less than 2 seconds, settling time of less than 10 seconds, and the steady-state error of less than 2%. For example, if the input is 0.2 rad (11 degress), then the pitch angle will not
• • • • exceed 0.22 rad, reaches 0.2 rad within 2 seconds, settles 2% of the steady-state within 10 seconds, and stays within 0.196 to 0.204 rad at the steady-state. Overshoot: Less than 10% Rise time: Less than 2 seconds Settling time: Less than 10 seconds Steady-state error: Less than 2% Transfer function and the state-space Before finding transfer function and the state-space model, let's plug in some numerical values to simplify the modeling equations (1) shown above. (2) These values are taken from the data from one of the Boeing's commercial aircraft. 1. Transfer function To find the transfer function of the above system, we need to take the Laplace transform of the above modeling equations (2). Recall from your control textbook, when finding a transfer function, zero initial conditions must be assumed. The Laplace transform of the above equations are shown below. After few steps of algebra, you should obtain the following transfer function. 2. State-space Knowing the fact that the modeling equations (2) are already in the state-variable form, we can rewrite them into the state-space model. Since our output is the pitch angle, the output equation is: Matlab representation and open-loop response Now, we are ready to observe the system characteristics using Matlab. First, let's obtain an open-loop system to a step input and determine which system characteristics need improvement. Let the input (delta e) be 0.2 rad (11 degrees). Create an new m-file and enter the following commands. de=0.2;
num=[1.151 0.1774]; den=[1 0.739 0.921 0]; step (de*num,den) Running this m-file in the Matlab command window should give you the following plot. From the plot, we see that the open-loop response does not satisfy the design criteria at all. In fact the open-loop response is unstable. If you noticed, the above m-file uses the numerical values from the transfer function. To use the state-space model, enter the following commands into a new m-file (instead of the one shown above) and run it in the command window. de=0.2; A=[-0.313 56.7 0; -0.0139 -0.426 0; 0 56.7 0]; B=[0.232; 0.0203; 0]; C=[0 0 1]; D=[0]; step(A,B*de,C,D) You should get the same response as the one shown above. Note: It is possible to convert from the state-space to transfer function, or vice versa using Matlab. To learn more about conversions, see Conversions Closed-loop transfer function To solve this problem, a feedback controller will be added to improve the system performance. The figure shown below is the block diagram of a typical unity feedback system. A controller needs to be designed so that the step response satisfies all design requirements. Four different methods to design a controller are listed at the bottom of this page. You may choose: PID, Root-locus, Frequency response, or State-space. Example: Modeling the Ball and Beam Experiment
M R d g • • L J r alpha theta gives can Problem Setup System Equations Matlab Representation and Open-Loop Response Problem Setup A ball is placed on a beam, see figure below, where it is allowed to roll with 1 degree of freedom along the length of the beam. A lever arm is attached to the beam at one end and a servo gear at the other. As the servo gear turns by an angle theta, the lever changes the angle of the beam by alpha. When the angle is changed from the vertical position, gravity causes the ball to roll along the beam. A controller will be designed for this system so that the ball's position can be manipulated. For this problem, we will assume that the ball rolls without slipping and friction between the beam and ball is negligible. The constants and variables for this example are defined as follows: mass of the ball 0.11 kg The design criteria for this problem are: radius of the ball 0.015 m lever arm offset 0.03 m gravitational acceleration 9.8 m/s^2 length of the beam 1.0 m ball's moment of inertia 9.99e-6 kgm^2 ball position coordinate beam angle coordinate servo gear angle Settling time less than 3 seconds Overshoot less than 5% System Equations The Lagrangian equation of motion for the ball is given by the following: Linearization of this equation about the beam angle, alpha = 0, us the following linear approximation of the system: The equation which relates the beam angle to the angle of the gear be approximated as linear by the equation below:
Substituting this into the previous equation, we get: 1. Transfer Function Taking the Laplace transform of the equation above, the following equation is found: NOTE: When taking the Laplace transform to find the transfer function initial conditions are assumed to be zero. Rearranging we find the transfer function from the gear angle (theta(s)) to the ball position (R(s)). It should be noted that the above plant transfer function is a double integrator. As such it is marginally stable and will provide a challenging control problem. 2. State-Space The linearized system equations can also be represented in state-space form. This can be done by selecting the ball's position (r) and velocity (rdot) as the state variables and the gear angle (theta) as the input. The state-space representation is shown below: However, for our state-space example we will be using a slightly different model. The same equation for the ball still applies but instead of controlling the position through the gear angle, theta, we will control alpha- doubledot. This is essentially controlling the torque of the beam. Below is the representation of this system: Note: For this system the gear and lever arm would not be used, instead a motor at the center of the beam will apply torque to the beam, to control the ball's position. Matlab Representation and Open-Loop Response 1. Transfer Function The transfer function found from the Laplace transform can be implemented in Matlab by inputting the numerator and denominator as vectors. To do this we must create an m-file and copy the following text into it: m R g L d J = = = = = = 0.111; 0.015; -9.8; 1.0; 0.03; 9.99e-6;
K = (m*g*d)/(L*(J/R^2+m)); %simplifies input num = [-K]; den = [1 0 0]; printsys(num,den) Your output should be: num/den = 0.21 ———- s^2 Now, we would like to observe the ball's response to a step input of 0.25 m. To do this you will need to add the following line to your m-file: step(0.25*num,den) NOTE: Matlab commands from the control system toolbox are highlighted in red. You should see the following plot showing the balls position as a function of time: From this plot it is clear that the system is unstable in open-loop causing the ball to roll right off the end of the beam. Therefore, some method of controlling the ball's position in this system is required. Three examples of controller design are listed below for the transfer function problem. You may select from PID, Root Locus, and Frequency Response. 2. State-Space The state-space equations can be represented in Matlab with the following commands (these equations are for the torque control model). m = 0.111; R g J = = = 0.015; -9.8; 9.99e-6; H = -m*g/(J/(R^2)+m); A=[0 1 0 0 0 0 H 0 0 0 0 1 0 0 0 0]; B=[0;0;0;1]; C=[1 0 0 0]; D=[0]; The step response to a 0.25m desired position can be viewed by running the command below: step(A,B*.25,C,D) Your output should look like the following: Like the plot for the transfer function this plot shows that the system is unstable and the ball will roll right off the end of the beam.
Therefore, we will require some method of controlling the ball's position in this system. The State-Space example below shows how to implement a controller for this type of system. If you are interested in knowing how to convert state-space representations to transfer function representations, and vice versa, please refer to the Conversions page. Autor: Zidarta [email protected] Perú 02/06/2008
| Página siguiente |