SOLUTION: a magnetic levitation system | My Assignment Tutor

Matlab and Simulink for Modeling and ControlRobert Babuska and Stefano StramigioliˇNovember1999DelftDelft University ofTechnologyControl LaboratoryFaculty of InformationTechnology and SystemsDelft University of TechnologyP.O. Box 5031, 2600 GA Delft,The Netherlands 11 IntroductionWith the help of two examples, a DC motor and a magnetic levitation system, the use of MATLAB andSimulink for modeling, analysis and control design is demonstrated. It is assumed that the reader alreadyhas basic knowledge of MATLAB and Simulink. The main focus is on the use of the Control System Toolboxfunctions. We recommend the reader to try the commands out directly in MATLAB while reading this text.The examples have been implemented by the authors and can be downloaded fromhttp://lcewww.et.tudelft.nl/˜et4092. The implementation is done in MATLAB version 5.3 and has also beentested in version 5.2.2 Modeling a DC MotorIn this example we will learn how to develop a linear model for a DC motor, how to analyze the modelunder MATLAB (poles and zeros, frequency response, time-domain response, etc.), how to design acontroller, and how to simulate the open-loop and closed-loop systems under SIMULINK.2.1 Physical SystemConsider a DC motor, whose electric circuit of the armature and the free body diagram of the rotor areshown in Figure 1.Figure 1: Schematic representation of the considered DC motor.The rotor and the shaft areassumed to be rigid.valuesConsiderforthe the following physical parameters:The input is the armature voltage V in Volts (driven by a voltage source). Measured variables are theangular velocity of the shaft ω in radians per second, and the shaft angle θ in radians.2.2 System EquationsThe motor torque, T, is related to the armature current, i, by a constant factor K: T = Ki.The back electromotive force (emf), Vb, is related to the angular velocity by:(1) . (2) moment of inertia of the rotordamping (friction) of the mechanical systemm2(back-)electromotive force constantK =0.01 Nm/Aelectric resistanceR = 1 Ωelectric inductanceL = 0.5 H VTJR L+Vb = Kωbω– + – 2From Figure 1 we can write the following equations based on the Newton’s law combined with theKirchhoff’s law:2, (7)and substitute it in (5) to obtain:. (8)This equation for the DCmotor is shown in theblock diagram inFigure 2.+ – (3).2.3 Transfer FunctionUsing the Laplace transform, equations (3) and (4) can be written as:(4)Js2θ(s) + bsθ(s) = KI(s),(5)LsI(s) + RI(s) = V (s) – Ksθ(s),where s denotes the Laplace operator. From (6) we can express I(s):(6) TorqueTs ()K Velocityωs1Js+b( ArmatureLoadBackemfVoltage AngleVs () θs1 s( KVsb()Ls+R 3Figure 2: A block diagram of the DC motor.From equation (8), the transfer function from the input voltage, V (s), to the output angle, θ, directlyfollows:. (9)From the block diagram in Figure 2, it is easy to see that the transfer function from the input voltage, V (s),to the angular velocity, ω, is:. (10)3 MATLAB RepresentationThe above transfer function can be entered into Matlab by defining the numerator and denominatorpolynomials, using the conventions of the MATLAB’s Control Toolbox. The coefficients of a polynomialin s are entered in a descending order of the powers of s.Example: The polynomial A = 3s3 + 2s + 10 is in MATLAB entered as: A = [3 0 2 10].Furthermore, we will make use of the function conv(A,B), which computes the product (convolution) of thepolynomials A and B. Open the M-file motor.m. It already contains the definition of the motor constants:J=0.01;b=0.1;K=0.01;R=1;L=0.5;The transfer function (9) can be entered in MATLAB in a number of different ways.1. As Ga(s) can be expressed as , we can enter these two transfer functions separately andcombine them in series:aux = tf(K,conv([L R],[J b]))Gv = feedback(aux,K);Ga = tf(1,[1 0])*Gv;Here, we made use of the function feedback to create a feedback connection of two transfer functionsand the multiplication operator *, which is overloaded by the LTI class of the Control System Toolboxsuch that is computes the product of two transfer functions.2. Instead of using convolution, the first of the above three commands can be replaced by the productoftwo transfer functions:aux = tf(K,[L R])*tf(1,[J b]);3. Another possibility (perhaps the most convenient one) is to define the transfer function in asymbolicway. First introduce a system representing the Laplace operator s (differentiator) and thenenter the transfer function as an algebraic expression:s = tf([1 0],1);4Gv = K/((L*s + R)*(J*s + b) + Kˆ2);Ga = Gv/s;It is convenient to label the inputs and outputs by the names of the physical variables they represent:Gv.InputName = ’Voltage’;Gv.OutputName = ’Velocity’;Ga.InputName = ’Voltage’; Ga.OutputName= ’Angle’;Now by calling motor from the workspace, we have both the velocity (Gv) and the position (Ga) transferfunctions defined in the workspace.3.1 Exercises1. Convert Gv and Ga into their respective state-space (function ss) and zero-pole-gain (function zpk)representations.2. What are the poles and zeros of the system? Is the system stable? Why?3. How can you use MATLAB to find out whether the system is observable and controllable?4 AnalysisThe Control System Toolbox offers a variety of functions that allow us to examine the system’scharacteristics.4.1 Time-Domain and Frequency ResponsesAs we may want plot the responses for the velocity and angle in one figure, it convenient to group the twotransfer functions into a single system with one input, the voltage, and two outputs, the velocity and theangle:G = [Gv; Ga];Another way is to first convert Ga into its state-space representation and then add one extra output beingequal to the second state (the velocity):G = ss(Ga);set(G,’c’,[0 1 0; 0 0 1],’d’,[0;0],’OutputName’,{’Velocity’;’Angle’});Note that this extension of the state-space model with an extra output has to be done in one set command inorder to keep the dimensions consistent.Now, we can plot the step, impulse and frequency responses of the motor model:figure(1); step(G);figure(2); impulse(G);figure(3); bode(G);You should get the plots given in Figure 3 and Figure 4.5Figure 3: Step and impulse response.4.2 Exercise1. Simulate and plot in MATLAB the time response of the velocity and of the angle for an input signalcos2πt, where t goes from 0 to 5 seconds.Figure 4: Bode diagram.5 Control DesignLet us design a PID feedback controller to control the velocity of the DC motor. Recall that the transferfunction of a PID controller is:6, (11)where u is the controller output (in our case the voltage V ), e = uc – y is the controller input (the controlerror), and Kp, Kd, Ki are the proportional, derivative and integral gains, respectively. A block diagram of theclosed-loop system is given in Figure 5.5.1 Proportional ControlFirst, try a simple proportional controller with some estimated gain, say, 100. To compute the closed-looptransfer function, use the feedback command. Add the following lines to your m-file:Kp = 100;Gc = feedback(Gv*Kp,1);Gc.InputName = ’Desired velocity’;Here Gc is the closed-loop transfer function. To see the step response of the closed-loop system, enter:figure(4); step(Gc,0:0.01:2);You should get the plot given in Figure 6:Figure 6: Closed-loop step response with a P controller.To eliminate the steady-state error, an integral action must be used. To reduce the overshoot, a derivativeaction can be employed. In the following section, a complete PID controller is designed.r DC Motoreω VelocityVVoltage+ –PIDFigure5:Closed-loopsystemwithaPIDcontroller.Time (sec.)Step Response0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 200.20.40.60.811.21.4From: Desired velocity75.2 PID ControlLet us try a PID controller. Edit your M-file so that it contains the following commands:Kp = 1;Ki = 0.8;Kd = 0.3;C = tf([Kd Kp Ki],[1 0]);rlocus(Ga*C); Kp =rlocfind(Ga*C); Gc =feedback(Ga*C*Kp,1);figure(9);step(Gc,0:0.01:5)The rlocus and rlocfind functions are used to select the overall gain of the PID controller, such that thecontroller is stable and has the desired location of the poles (within the defined ratio among the Kp, Ki and Kdconstants). If the design is not satisfactory, this ratio can be changed, of course. We should obtain a plotsimilar to the one in Figure 7:5.3 Exercise1. Use the root locus and the Nyquist criterion to find out for what value of the gain Kp the proportionalcontroller for the angle Ga(s) becomes unstable.Step ResponseFigure 7: Closed-loop step response with a PID controller.6 SIMULINK ModelThe block diagram from Figure 2 can be directly implemented in SIMULINK, as shown in the figure Figure8:Time (sec.)0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 500.20.40.60.811.21.48Figure 8: SIMULINK block diagram of the DC motor.Set the simulation parameters and run the simulation to see the step response. Compare with the response inFigure 3. Save the file under a new name and remove the position integrator along with the ‘Graph’ and ‘ToWorkspace’ blocks. Group the block describing the DC motor into a single block and add a PID controlleraccording to Figure 5. The corresponding SIMULINK diagram is given in Figure 9. Experiment with thecontroller. Compare the responses with those obtained previously in MATLAB.7 Obtaining MATLAB Representation from a SIMULINK ModelFrom a SIMULINK diagram, a MATLAB representation (state space, transfer function, etc.) can be obtained.The ‘Inport’ and ‘Outport’ blocks must be added to the SIMULINK diagram, as shown in Figure 10.Then we can use the linmod command to obtain a state-space representation of the diagram:[A,B,C,D] = linmod(’filename’);where filename is the name of the SIMULINK file.Figure 9: SIMULINK block diagram of the DC motor with a PID controller.Figure 10: SIMULINK block diagram of the DC motor with ‘Inport’ and ‘Outport’ blocks.omegaTo WorkspacetClock To Workspaceangularspeed 1 s Step Input Ls+RK(s) Armature Js+b1 Load+ –K thetaTo WorkspaceangleDC motoromegaTo WorkspaceangularspeedStep Input + – Sum PID ClocktTo Workspace 1 s Ls+RK(s) Armature Js+b1 Load+ –K1Outport1Inport97.1 Exercise1. Convert the four matrices A, B C and D into a corresponding state space LTI object. Convert this oneinto a transfer function and compare the result with the transfer function entered previously inMATLAB.8 Linearization of Nonlinear Simulink ModelIn this section, we will use an example of a highly nonlinear system to show how to linearize a nonlinearSimulink model and extract the linearized model to MATLAB for control design purposes.8.1 Magnetic Levitation SystemMagnetic levitation as a friction-less support for high-speed trains, in bearings of low-energy motors, etc. Itconsists of an electromagnet which is attracted to an object made of a magnetic material (such as a rail).The control goal is to keep the air gap between this material and the electromagnet constant by controllingthe current in the coil. A schematic drawing is given in Figure 11.Figure 11: Schematic drawing of the magnetic levitation system.The position and the motion of the object in the magnetic field are dependent on the forces that act on it.These forces are: (i) the gravitational force, (ii) the electromagnetic force, and (iii) a disturbance force. Thedynamic equation of the system is derived from the basic law F = ma,, (12)where Fgrav = mg, Fdist is an unknown disturbance, and the electromagnetic force is. (13)In our example, we use Kmag = 17.8 µH and m = 8 kg.8.2 Nonlinear Simulink ModelThe nonlinear equation (12) is implemented in a Simulink model given in Figure 12.RailAirgap zCurrent iFRFgrav FdistComputerD-AA-D10Figure 12: Nonlinear Simulink model (bearing.mdl) of the magnetic levitation system.8.3 LinearizationLet us linearize the nonlinear model around an operating point y0 = 2 mm. There are two possibilities tolinearize a nonlinear model:• Analytically: by hand or using symbolic maths software such as Mathematica, Maple or the SymbolicToolbox of MATLAB.• Numerically by applying the trim and linmod functions of MATLAB.The second possibility will be explored here (you can do the first one as an exercise). Let us use thefollowing script (lin.m): params;% a script with definition of system’sparametersfile = ’bearing’;% nonlinear Simulink model to be linearizedu0 = [10; 0];% initial input guess [input; disturbance]y0 = 0.002;% initial output guessx0 = [y0 0]’;% initial state guess [x0,u0]=trim(file,x0,u0,y0,[],[2],[]);[A,B,C,D] = linmod(file,x0,u0); sys =ss(A,B,C,D);% make an LTI object The trim function numerically searches for an equilibrium of thenonlinear system. A reasonable initial guess (x0, u0 and y0)must be provided. The additional parameters of this functionare indices of the inputs, states and outputs that are not free tovary during the search. A typical example of such a variable isthe state variable corresponding to the operating point.The linmod function extracts the matrices of a linear model obtained by numerical linearization at theequilibrium. Once this model is available, it can be used for analysis or control design.y’’ y’1 sCheck limits Air gap 1 1 s u2 Kmag/m1/m 1/u^2 g 2Disturbanceforce1Current118.4 Exercise1. Choose another operating point and extract a linear model at that point. Compare to the modelobtained above in terms their gains, poles and zeros.9 Concluding RemarksThe authors hope that this text has been useful and would appreciate receiving feedback from you. Let usknow if you have found any errors and omissions or if you have suggestions for improvements. Send thempreferable by e-mail to: R.Babuska@ITS.TUDelft.NL.