About ENGM 541, ENGM 670-X5, and MECE 758-X5
Modeling and Simulation of Engineering Systems
Winter 2011
Instructor: M.G. Lipsett
This FAQ is maintained by the course instructor, and is located at http://www.ualberta.ca/~mlipsett/ENGM541/FAQ.htm
It was re-organized on January 7, 2010.
IMPORTANT LINKS
The course website is located at http://www.ualberta.ca/~mlipsett/ENGM541/ENGM541.htm
Course Administration Questions
Questions Related to This Year’s Assignments
FAQ Questions from Previous Years
If a student asks a question by email, then the answer gets posted here
on the FAQ, so that all students have fair access to any additional information
that may be provided. Other information that will be of interest to the whole
class gets posted here as well.
13) How long does the
project report have to be?
The report should be as long as it needs to be. Typical reports are in the 20 – 30 page range (plus references and any appendices), but that’s just a guideline. Page length will vary depending on the topic, because it will depend on how long it takes to explain your model and your simulation results. Remember that your report must include any source code on an accompanying CD (labeled). A good concise report is better than a long mediocre report.
If you are developing your own governing equations, then a simple simulation for a few trial sets of parameters and input conditions would be sufficient to demonstrate that you have attempted to verify the simulation. A literature review will have to be well annotated and of course fully referenced. A simulation based on an existing set of governing equations will have to have good verification and preferably validation for the application of interest, and a number of scenarios to show the behaviour of the model and the sensitivity of the outputs to different parameters and input conditions. A probabilistic simulation will require justification for the distributions used. (Apr 6/11)
12) What will be on the 2011 final?
The final exam is cumulative, meaning that all material presented in the course may be tested, including lectures, labs, related text materials, assignments, & resources posted on the course website. The final will be held on Wednesday April 13, 2011 during the regular lecture time, 5 pm to 7:30 pm, in the regular lecture room ETLE 2-001. There may be questions about numerical procedures and related MATLAB methods; but there will not be any questions that require access to a computer to solve. Questions will be appropriate for calculations by hand and by calculator. The exam will be two hours-and-a-half hours long, open book and open notes, calculators allowed (but communication features must be turned off). No computers. Sample Final (with solutions) is posted in the Examples directory on the course web-site; it is fairly representative of the number and type of questions that might be asked, although the specific number of questions and format of questions may be different. Make sure that you show your work, so that understanding is clear about how to formulate problems with appropriate choices of variables that satisfy admissibility requirements, what constitutive relationships need to be used, and solution methods for particular types of problems. Marks are given for demonstrating understanding as well as for the correct answer, so make sure that you show all the steps in your work. Graduate students (ENGM 670 and MECE 758) can expect some more challenging exam questions than the ENGM 541 exam questions.
11) Can I submit the project electronically?
Unfortunately, no, unless you are granted special permission in writing. Once the projects are graded, they will be available outside the instructor’s office. In some cases you may be asked if a copy of the project can be kept for later use in the course. If any project material is used in the future, then credit will be given for the source of the material. (Mar 31/10)
10) How do I submit
the code for my project?
If you have developed a simulation, then make sure that you include the source code. Marks will be deducted if code is not provided. A CD is best. I prefer that you not send it as an email attachment. A hard copy as an appendix is acceptable, provided that all the required files are included. If any project material is used in the future, then credit will be given for the source of the material.
9) I am doing a discrete-event model for my project. I do not have access to SimEvents, and my system looks kind of complicated to do in Excel. There are some other packages available for doing DES, so can I use Extend, Arena, Ideas, Stella, or AweSim for my project?
Yes, with written agreement from the instructor you can submit a project done using a different modeling package. It is doubly important that your project gives a complete description of the model formulation and assumptions, and that you provide source code for your project (CD format, not email). (Mar 4/11)
8) When is the project update due? What is the format?
The project update is due on March 23rd 2011. There is no set format for the update, other than it has to be at maximum of five pages long, and submitted as a hard copy (email is not sufficient). The update is not marked. It can be a draft of your report, or an outline of the formulation of your governing equations, or whatever you want. If there are issues with your scope, please identify them and discuss how you will modify your plan to achieve your objective )or modify your scope). The intent is that you show some progress, NOT dazzle the instructor with a glossy report. You can give it to the TA in the lab, or drop it in the assignment box on the fourth floor of the MECE building.
7) Is there an
example or a template for the project?
The instructor has a number of projects from previous years that students have given permission to use. There is no specific format for the project. You are free to use whatever format works for you provided that you have the following sections in your report: introduction, project objective, definition of the system or technique of interest, methods used (literature search, analysis methods, programming), results, discussion, interpretation, and references. The project will be marked for presentation quality as well as technical quality.
6) What topic should
I propose for the project, and what is expected for the project?
You are free to propose any topic in modeling and simulation, subject to approval by the instructor. So choose something interesting. The instructor has a list of ideas if you’d like suggestions. (They are in the Projects folder on E-Class, and there’s a link on the course website.)
There are three general types of project that you might consider:
a) Define a problem in engineering analysis for a (small) physical or technological system, generate the governing equations, define a set of trial parameters based on well-described assumptions, generate simulation results, and discuss the results and their physical interpretation.
b) Research a topic in systems modeling methods and generate a summary report describing the method and give examples of its applicability. Candidate topics include bond graph theory, linearisation, system identification, sources of errors in modeling and simulation, nonlinear eigenvalue systems, random process modeling, machine learning, modeling dynamic financial systems, causality (directed graphs), dynamic programming, animation techniques, human-machine interaction using simulators (such as flight simulators, nuclear reactor simulators, surgical training simulators), chaotic systems, using modeling and simulation for control, reliability, optimization, etc.
c) Take a known problem in systems modeling and develop a Matlab/Simulink simulation of the system, with documented code, and a report comparing your simulation results to some benchmark (published experimental data, analytical results, etc.) Systems should be somewhat complex, involve physics that have not been considered in the course so far (e.g. optics, reactive processes, etc.) or mixed systems such as microelectromechanical systems (MEMS) or robots. Developing a solution in Matlab/Simulink should entail an interesting simulation technique such as animation.
The instructor may ask you to modify your proposal if the scope is too broad (or too narrow). Projects can not be submitted on work that has already been done for previous courses or for thesis work. Projects may, however, contribute to thesis work that has yet to be done. All projects require a report from each student and any source code that was developed (such as m-files). Each report must include an introduction, project objective, definition of the system or technique of interest, methods used (literature search, analysis methods, programming), results, discussion, interpretation, and references. Source material must be referenced, including web sources. Some Matlab/Simulink code is available for open source use on-line; and this access to existing code is a huge advantage for the Matlab community. Use of such material is acceptable but it must be referenced, so that your contribution is clear.
The instructor may wish to use excerpts from reports or sample code in future offerings of the course, in which case the instructor will ask for permission. If permission is granted, then the source of any material used will be acknowledged.
5) What course
material will the midterm cover, and what’s the
format?
The midterm will be held on Wednesday
March 2, 2011 during the regular lecture time, 5 pm to 7 pm, in the regular
lecture room ETLE 2-001. The midterm will test formulating governing
equations and methods of solution. There may be questions about numerical
procedures and related Matlab methods; but there will
not be any questions that require access to a computer to solve. All material
presented up to and including Feb 17, 2010 (lectures, laboratory handout
information, posted material on course web site & e-class) may be tested.
Questions will be appropriate for calculations by hand and by calculator. The
midterm will be two hours long, open book and open notes, calculators allowed
(but communication features must be turned off). No computers. A Sample
Midterm (with solutions) and Another (older) Sample Midterm (with solutions) are
posted (with solutions at the back) in the
Examples directory on the course web-site & on E-Class. They are
representative of the number and type of questions that might be asked,
although the specific number of questions and format of questions may be
different (in previous years eigenvalue systems were
covered earlier in the course; you won’t be responsible for eigenvalue
problems on this midterm). Make
sure that you show your work, so that understanding is clear about how to
formulate problems with appropriate choices of variables that satisfy
admissibility requirements, what constitutive relationships need to be used,
and solution methods for particular types of problems. Marks are given for
demonstrating understanding as well as for the correct answer, so make sure
that you show all the steps in your work. Graduate students (ENGM 670 and MECE
758) can expect some more challenging exam questions than the ENGM 541 exam
questions.
4) Is attendance at
the lab sessions mandatory?
For ENGM 541 students, the labs will be marked (mostly for effort). Attendance is not taken at the lab sessions; but it is highly recommended that you attend. Material that is discussed during labs may likely appear on the exams. Material presented in labs will be posted on the course website.
3) Is the textbook
mandatory?
The course textbook is recommended but not mandatory. We will be making extensive use of the text in the beginning and the middle part of the course, particularly concepts of how to build models of lumped-parameter systems for different types of physical systems, and the dynamics of such systems. Some problems will come from the text. The book has been ordered for purchase at the bookstore. Crandall’s book is also worth buying as a resource (it is available used on-line), and some of the other listed resource books are also very worthwhile.
2) Are the notes for
the lectures available?
Hard copies of the lecture notes that are posted on the course website will
be provided so that students can do their own annotations. The lecture notes are available for
downloading from the course web-site & e-class. Annotated slides are posted
and audio files may be posted as well.
1) How do I print the lecture slides?
One of the reasons for posting the slides on the website early is so that students can look at the material in advance of the class, or make their own copies to annotate. Not all of the lecture information is printed on the slides - there is blank space provided for additional notes. Please note that these slides are copyrighted, and can not be distributed or modified without the copyright holder’s permission. The slides that are posted so far for the lectures are in pdf format with two slides per page. (Last year the slides were also available in other formats but students preferred this format.) On a PC, you can view them directly in Internet Explorer with the Adobe plug-in, and then print from there. You can save a copy on your own computer and use Acrobat Reader to view and print the slides (respecting copyright, of course). If you have trouble printing, first make sure that you have the latest update of Reader. If all else fails, you are welcome to stop by my office to borrow a set of hard copies of lectures to date, which you may photocopy once for personal use, and then return my set. Please let me know if you have continuing trouble printing the slides, or if the URL is a broken link. (Dec 13/08)
ASSIGNMENT #5:
There is a typo on
Problem #2. Units of L are demand/unit/day^2, or demand x
days^2 /unit. This will give you correct frequency units of rad/day.
(Other parameter units are correct.)
I am working on Assgt 5 problem 2 and I’m not sure how to find the damped
frequency and plot the response.
This system is a linear second-order system, and so there are analytical methods in the text for finding the system response. There are two basic choices for generating the response spectrum.
1) If you have the state-space representation of the system in matrix form (or the transfer function), then Matlab can generate the natural frequency and damping coefficient with the command damp (from which you can find the damped natural frequency) and you can produce a bode plot of the ‘free response’ of the system using the command bode. There is a worked example in the posted solutions from previous years that shows this, as well as examples in the text (based on mechanical systems, but the idea is the same). To produce an analytical solution such as a bode plot, you will require some initial condition. You can use initial values of 0 for the variable of number of units and 1 for the variable of rate of units transported in the governing equations.
2)
There is a way to generate the general
frequency response function without having to have initial conditions; but to
be honest it requires a fair amount of work. You can select an appropriate Nyquist frequency, and thus determine the time step delta
t. You also have to choose a length of sequence that makes sense for the
conversion to the frequency spectrum; for the FFT a power of 2 is needed – so
512 or 1024 points is usually a good number. Then “drive” the system with an
input that excites all of the possible frequencies of the system (Gaussian
noise), without putting in an excitation frequency that can be aliased (this
would generally mean putting the input through a low-pass filter). Then, once
that’s all done, the input sequence is used to find the time series response in
a simulation using an ODE solver. The output from the solver is the output of
the system that is now ready to be analysed, and then
plotted with respect to frequency, using the FFT or PSD. As you can tell, this
requires a lot more work; but if the system is nonlinear (and does not have a
transfer function or a linear state-space formulation), that’s what we have to
do, unless we linearise around a normal operating
point. (Mar 25/11)
I have a question
regarding the vibration analysis of systems where the motion is in a direction
parallel to the force of gravity. When analyzing these systems I have
noticed that in the textbook, as well as in the notes, we typically ignore the
gravitational force term ‘mg’. I understand this could be a reasonable
assumption when this force is small compared to the forcing function, f(t), or the spring/damper forces, but is this something
that is commonly done in practice? I certainly realize it simplifies the
problem when we neglect the gravitational force, but it just seems odd that
when analyzing these systems we wouldn’t take it into account….
We always do our loop equations in terms of a change in the loop variable. We could start from an undeflected spring, and include a gravity force, thereby including all forces that act on the system. But if the force is constant and the node (mass) is acted on by an equal and opposite force, then we can ignore that bias, thereby allowing us to avoid having to write an expression for any force that is already compensated in the system. If we have a mass hanging off a spring, and the spring already has a reference displacement delta_0 that equals the weight of the mass, then that cancels out in the problem. A Mass-spring system under the same input forcing function will have the same oscillatory motion in the vertical direction whether it is on earth or on the Moon; the only difference will be the reference displacement of the spring due to the weight.
If the system is nonlinear, then we usually will consider all forces, because the normal operating point is affected. Hope this helps. (Mar 25/11)
Is there any extra
reading material you could recommend for further learning on discrete-event
simulations?
An excellent general guide to modeling and simulation is the Handbook of Systems Engineering and Management, by Sage & Rouse, which is available through the U of A library electronically. Also, Simulation Modeling and Analysis (3rd Edition) by A.M. Law is probably the best generally text on discrete-event modeling in the field. Reliability modeling is covered, but not in huge detail; Kuo & Zuo’s book is good on that special class of modeling. If you are interested in systems thinking, then you should start with The Fifth Discipline by Senge; but then move to the works of Russell Ackoff – they are terrific. And for insights into how quantitative analysis translates into good morale and high quality, read W. Edwards Deming. (Mar. 24/11)
ASSIGNMENT #4:
I'm a little confused
by the diagram for problem 4 of assignment 4.
Are we supposed to treat this as a buckling problem where the
displacements are in the horizontal direction as shown in the diagram? Or would the masses just move up and down in
the vertical direction?
The masses oscillate with a horizontal motion. The vertical columns act as cantilever beams that produce a horizontal spring force when there is relative horizontal displacement between floors. I hope this makes the problem a little clearer. (Mar 12/11)
For part c) on
problem # 7 I don't see how the expression is any different than found in part
b). As drawn in figure 2 beta looks like
it equals theta and so they are the same angles, resulting in the same
expression for the damage rate, unless I am missing something here?
Problem 7 part c is asking for the damage rate as a function of theta around the part of the circumference between - beta < theta <= + beta. Theta does not equal beta. (Mar 12/11)
I was working on
Assignment 4, question 1, and I have a question. I'm pretty sure I've derived
the second order differential equations, however I do
not know how to derive the transfer functions... In the text they use Laplace
transforms, but do not show their work as to how they get to the Laplace
transform. Also, does the block diagram need to be Laplace variables, as shown
in the examples of the text? I'm guessing not, since lower case omega_L is the output.... Finally, just for clarification
on block diagrams, since I've never done a multiple differential equation one,
will it require three loops for the three equations? I was hoping to try all
parts a,b,and c to make sure
I understand the material. Thanks.
The idea of generating equations based on Laplace transforms is pretty straightforward. Recall that in our general framework, we generate admissibility equations based on one set of variables and then use the constitutive relationships to express the governing equations in terms of the variable type in which we want to express the governing differential equations. To generate Laplacian expressions, we construct admissibility equations in terms of one variable type; then -- rather than substituting for the complementary variable using a constitutive relationship that is in terms of the complementary variable (or an integral expression or differential expression of the complementary variable) as we ordinarily do -- we substitute the appropriate Laplace transform for each element, and for any input function (which we ordinarily think of as the constitutive relationship for a source). There are some examples in the text of simple systems expressed with Laplace transforms, for which the differential equations are easy to check. There is also at least one worked example posted that has multiple loops in a Simulink block diagram, which is the same idea. There will be as many differential equations as there are unique loops in the system, and each differential equation will be the order determined by the elements in the loop (first-order, second-order, etc.) I hope this helps. (Mar 13/11)
ASSIGNMENT #3:
Hi Dr. Lipsett, I
have some questions concerning the robot arm stiff example in lecture 5,
p.16-17. I'm having a difficult time
getting the governing equations that you found. I also am having a hard time
connecting the figure on slide 31 vs the figure on
slide 32, and seeing where they match up. What are the components shown in both
of the figures? I think I got the first equation on slide 33 in the correct
way. We are doing a sum of moments correct? About what point though, and where
do we get the variable l from? Also, is the moment of inertia term considered
opposite in direction in comparison to the input torque on the very most left
box? (I just used the variables you had in the blocks which included the length
for the moment sum.) Also, in the second equation on slide 33, I don't get the q_double_dot in the formula. Could you clarify where that
comes from, please? Same with the y_double_dot
in the third equation. Thanks.
There are three rigid bodies, but only two of them are always part of the arm. One of them m_e is the workpiece that the force sensor on the end of the arm may contact. l is not a variable; it is the length of the arm. We are using the small angles assumption for the moment balance, so the diagrams are representing small motions of the masses. q is the displacement with respect to an inertial frame for the connection between the end of the arm and the end effector. q double dot is an acceleration term dictated by loop admissibility by differentiating the original loop equation twice. Same for y double dot. I suggest that you give a listen the segment of the audio file for this lecture where I set up the problem. Hope this helps. (Feb 19/11)
On Lecture 5 slide 39
I noticed when I was posting the lecture that I had thrown in an extra negative
sign for the negative effective torque term of equation 2. It is correct on the
annotated slide deck. Sorry for the slip of the pen. (Feb. 19/11)
I was wondering how
much of the extra grad level material will be on the midterm. I’m completely
lost in the chaos… Also, I have tried
simulating the chaotic system, but I don’t get exactly the plot that is shown.
It’s easy to get lost in phase space, especially in a chaotic system…
In many engineering systems, we design to avoid interesting behavior, and change variables slowly, which is why equilibrium analysis is often enough. The whole point of chaotic systems is to recognize that they are nonlinear dynamic systems that have particularly non-intuitive behaviour – but that modeling and simulation can help us to understand the nature of these systems. For local regions, we can often get sufficient insight by linearising the system to determine the locally linear behavior, which we will explore in our next set of lectures on eigenvalues (with extra info on singular value decomposition) and spectral analysis techniques.
The intent of the extra material is to expose graduate students to more sophisticated problems and analysis techniques than the undergrads get. Because we don’t cover it in the depth that the regular material gets, it will not be a primary focus for examination material.
As for the phase portrait of the pendulum system with the two magnets, you will expect to see a funny oscillating behaviour, but no spiraling down to an attractor, because the system that we are simulating does not have a dissipating term that causes energy to leave the system. The specific shape of the oscillation will depend on the parameters that you choose for the system. (Feb 18/11)
Regarding problem 4,
can you please clarify what you mean by "Run the simulation using an Euler
integration scheme"? When I check out Euler in the notes, it is a
numerical method for solving initial value problems (Lecture 5 slide 15).
However question 4b doesn't have any initial values. The Euler method problems
I found in previous assignments also included initial values. Are we to assume
0, 0, 0?
You are correct. You need to choose some initial values to start the simulation. Different initial values will lead to different trajectories in each solution. (Feb 17/11)
I know the assignment
is due on Friday, but you said something in class about being able to hand it
in late.
This assignment is shorter to reflect the fact that not everyone will be on-campus during Reading Week; but you may submit the assignment into the marking box on the 4th floor of the Mechanical Engineering Building up to 5 p.m. on Wed. Feb.23, 2011.
How do I use an Euler
integration scheme in Matlab? I've generated
results using ode45, but when I replace that function with ode1 (Euler), I get
an error - undefined function.
ODE45 is an R-K routine. You don’t use it for Euler integration. For Euler integration, you do your own, fixed-step integrations, evaluating f(x,t) at time t, and updating the estimate of x at time t + delta_t, then moving forward a step, and doing it again until you get to the end time t_f. This is easily done using a for loop. If you google “euler integration matlab”, there are a pile of examples out there. But make sure you create your own version! (Feb 11/11)
ASSIGNMENT #2:
In Wednesday's class
you did a least squares example to solve for undetermined coefficients (slide
47, example of least squares 2). The equation resulting from differentiating
with respect to c1 was non-linear in both c1 and c2 yet the next slide arranges
the (hypothetical) equations in a matrix. Is the resulting system supposed to
be solved iteratively, or will the resulting equations always
reduce to a linear system?
In methods of undetermined coefficients, we create expressions in terms of the coefficients ci and then solve for these coefficients. Looking at the equations again, you'll see that the expressions are linear in terms of the coefficients (although the expressions in the elements of the A matrix are complicated, they are in terms of constants), which is why we can put them into matrix form to solve for the vector of ci terms. (Feb 5/11)
ASSIGNMENT #1:
Note that Problem 4b is
problem 1.10 from the course text (not 1.11). This has been reposted on the
assignment. If anyone does 1.11 it will also be marked (but 1.10 is easier).
I am having trouble
coming up with the node and loop variables for question 2a. I am a bit confused with this problem as when
finding the node and loop equations I end up with 3 node equations and 2 loop
equations for only 4 unknowns.
You are on the right track. Really. When you look at the node equations that you have generated, ask yourself whether the third one actually gives you any additional information not already provided by the first two admissibility equations. (Here’s a hint: the answer is no.) You will recall from class that once we have satisfied node admissibility (which gives us expressions for the two node variables that we will eliminate in the final set of equations in terms of the variables that we will use), then we turn to loop admissibility to provide us with a sufficient number of equations to generate a solution. Think about how many node variables you have left over after satisfying node admissibility, and also think about how many unique loops to requires to cover all of the elements in the system. I think that you will find that you don’t actually have a problem, just an extra equation that you don’t really need, and which doesn’t add any extra information to contribute to the solution. (Jan. 20, 2011)
EQUILIBRIUM
SYSTEMS
For the direct
formulation for governing equations, which is the best variable type to use?
The choice of variables often depends on how you want to describe the state of the system. If you have a choice, loop variables are usually easier, because it is simpler to satisfy the admissibility requirements. Provided that you define the loop variables at the nodes, the loop variables are unique and single-valued, guaranteeing loop admissibility. (But you should always state that in your solutions.)
In Lecture 1, on
slide 30 (pg 15), we identified some loop and node variables for four types of
systems. Can you provide any analogous
examples that relate to systems with an ENG M focus? I think financial transactions or production
rates might be examples, but further insight would be appreciated.
Technological systems (and business processes) generally deal with flows and transformations of money, material goods, information, energy, and people. The concept of loop and node variables related to energy transfer within physical systems extends to financial systems by using flow of money as a node variable and utility (how much value the element adds to the system) as the loop variable. As we will see later in the course, utility for part of a system is difficult to assess, and it depends on the business context. Many models also do not use explicit loop and node models, but rather use one variable type only, with some rules for calculating how the system behaves according to that variable type. Most financial system models, inventory models, scheduling models, etc. take this approach. Using discrete-event systems allows those rules for how elements of the system behave, because we only care about what happens (and when) the system changes. (Jan 15/11)
Why is co-energy U*
negative when using a voltage source?
A source puts energy into the system, which means that energy is coming out of it. In other words, a source is a driver, not a load, within the system. By our convention, if energy is being stored within the element, then the term for that element in the extremum formulation has a positive value. Since a voltage source is releasing energy into the system, its term has a negative value.
What do I do when I
am having trouble nondimensionalizing the equations
and I can’t get rid of one term?
There are often different ways to nondimensionalise a set of equations. Depending on how the loop and node variables are defined, and how the nondimensional variables are defined.
In this problem, we have governing equations in terms of temperature, and so we create nondimensional variables by dividing by some temperature (not a temperature variable though, because that would produce a nonlinear set of equations!). We are provided with given information about the (constant) forcing temperatures in the problem: T0 and T4. We could use either one of these constant parameters to nondimensionalise, or we could even use an expression in terms of both of them, such as (T0 + T4 ) or (T0 - T4).
We can also express our loop variables in different ways. For example, you could use the temperatures T1, T2, and T3, themselves as loop variables, but you could use another set of temperature expressions as well, such as (T1 - T∞), (T2 - T∞), and (T3 - T∞), where T∞ is some reference temperature (such as absolute zero). But using a new reference parameter T∞ would not be very convenient. Better to use one of the parameters that is given in the problem (T0 or T4). If you do that, then note that all of the loop variables need to be expressed with respect to the reference loop parameter. The result is that any term that includes only the constant parameter used as a reference will cancel out. By using this trick, you should be able to eliminate that annoying term and produce the set of nondimensional equations. If you look one of the example problems posted on the course website (a direct solution to the pipe network problem we solved by an extremum method in Lecture #3), you will see how this is done. (Jan 15/10)
What are the normal
operating conditions for linearising?
Use steady-state conditions (derivatives are all equal to zero), and then find the values of the loop and node variables (e.g. voltage and current) through the nonlinear element under those conditions. From the normal operating point, use incremental variables to find the linearised form of the constitutive relationship for the nonlinear element.
Any governing equation that includes a term related to this element has to be rewritten, using the linearised constitutive relationship for the element in place of the original nonlinear constitutive relationship.
We have two options for doing this. The preferred approach is to rewrite the governing equations in terms of the incremental variables, including the substitution for the linearised constitutive relationship. This yields the linearised state equations about the normal operating point.
A more complicated approach is to find the linearised form of the constitutive relationship, and then substitute back to express the linearised constitutive relationship in terms of the regular variables (ec, i2), not the incremental variables. This substitution will also yield the linearised form of the governing equations, because there is only one nonlinearity in the system. Here is the explanation in pdf form for the second option, including the linearisation equations.
Take a look at lab #5 for some hints on the linearisation procedure, or look at the course text. (Feb 22, 2010)
Why don’t we do all the
possible directions of motion when we do the admissibility equations for a
mechanical system?
When we do our loop equations, generally we have to consider all the possible directions of motion (three translational and three rotational for each body); but for specific problems we can make life simpler by neglecting to do loop equations for which no admissibility equations can be written; that is, there is no displacement for which there is an element providing energy or receiving energy from the system.
On 2010 Assignment #4
Problem #4, I’m having trouble getting the transfer function form from the
system state-space form. When I use ss2tf I get an error message saying that there
are not enough input arguments. What can I do?
(UPDATED ANSWER) This system has two transfer functions, because there are two outputs of the system: x1 and x2. There is only one input, and so the transfer functions are X1(s)/F(s) and X2(s)/F(s). The command ss2tf produces only one transfer function, and if the system can potentially have more than one transfer function, then ss2tf requires an additional argument. If you type help ss2tf at the command line, then you’ll get information on the syntax, which is ss2tf(A,B,C,D,iu) where A,B,C, and D are the state space matrices, and iu is a scalar number denoting which input is to be used. Here, there is only one input but you still to include the argument. If you don’t assign ss2tf to a vector, then the command will return an output that is a matrix of 2-element column vectors representing the numerators of the two transfer functions (because in this case there are two outputs).
A simpler approach is to use the ss command (which takes A,B,C, and D as arguments), assign that to a vector representing the system, and then use that vector as the argument in the command tf to yield the first transfer function form, and then use that result to produce the response of the system using the damp command. Although this system has two degrees of freedom, it may not necessarily have two damped frequencies.
If you assign a pair of vectors for the command, e.g., [Y,X]=ss2tf(A,B,C,D,1) then you will get in Y the numerators of the transfer functions (two in this case), and X is the vector of coefficients of the denominator for each of the transfer functions. X is only a row vector because there’s only one input. To apply the damp command correctly for the pair of matrices X and Y, check the Matlab help, because you’ll only use one row of the matrix Y to produce all of the eigenvalues of the system and the associated natural frequencies and damping ratios. It’s tricky to apply this correctly to get the overall system response. I recommend the approach outlined in the preceding paragraph. (March 13/10)
I'm not sure of how
to go about finding the transformation matrix that diagonalises
[A]. The system is simple enough that I can get the characteristic polynomial
to find the eigenvalues and corresponding
eigenvectors... can I just do that?
The question is asking that you use the diagonalisation method of Jacobi rotations. It's a good insight on your part that by using an alternative method it would be possible to find the eigenvalues and corresponding eigenvectors, and then construct the diagonalised representation of the system, but that's not what the question is asking. Hope this makes it clear. (Jan. 30/09)
(Supplemental) I gave
that a try but I'm not sure how to deal with the 3x3 [A] matrix and the 2x2
transformation matrix. The way I'm doing it, I just add more non-zero terms on
the off-diagonal and some of my diagonal terms become 0.
The two by two transformation matrix is a "submatrix" that is used in the "big" transformation matrix that has the same dimensions as the H matrix (which in the assignment problem is 3x3). If you're trying to zero out elements h_pq and h_qp in the current version of the H matrix, The only nonzero elements of the big transformation matrix will be T_pq, T_qp, T_pp and T_qq. Take another look at the example spreadsheet that I went through at the end of lecture 4; it shows how the 2x2 T submatrix gets applied. It should appear with its diagonal on the diagonal of the big matrix. (Feb. 1/09)
I started the
homework for the first assignment, and I think this might be the first email of
many. I’m not sure how to get started. Will you be in your office today?
For general set-up of the solutions to equilibrium problems, You should take a look at the example problems found in the Examples directory on the course web site, including the solutions for the first problem set in 2008, which are posted. Follow all of the steps to be sure that you have the correct equations for admissibility. You should not necessarily expect to get the same set of governing equations by both the direct and extremum methods (but if you solve the governing equations the solution vector will of course be the same).
In 2008 assignment 5-
Q5, why are the directions of Torque & speed for the turbine in the
solution different with what has been shown in the book for the same system? (P
254 of the book)
Torque is delivered by the turbine to the mechanical system (shaft, friction, and mechanical load of the generator). In the solution, it is shown as being the driver into the mechanical system. For that reason, the torque Tt and the direction of motion omega are in the same direction (being the driver). I don’t have my copy of the text with me, and so I can only guess that in the figure it is showing the turbine as having flow of water going in and producing (driven) torque out of the ideal turbine which would be in the opposite direction to omega (and which then becomes the driving torque input into the mechanical subsystem). In any diagram it is important to represent what is providing power to a subsystem and how the power is being consumed. Hope this helps. (Apr 9/08)
For solving a
differential equation, is that also valid to transform it to first-order
equations (i.e. equation set) and proceed? This is dealing with 2008 Assignment
question 3.
While it is certainly possible to solve a Taylor Series problem by breaking a differential equation of higher order into set of standard first-order is possible, it offers few advantages, because in each case you would still need to use successive differentiation. Note that Picard’s method can be extended to differential equations of higher order, but it is conceptually much easier to generate a solution by putting a system involving higher-order differential equations into first-order form. (Apr 9/08)
In 2008 assignment 5
question 4, I just wanted to check the definition of utility. Is the following
math relation correct: Utility=(income-Cost)/Cost ?
Generally, benefit (also called profit or income) is the integral of cost with respect to delta utility. In question four, you can treat the cost as a “source,” and so for our purposes in this problem we can say that delta utility is profit/cost. Profit is also called income, and income is revenue minus costs. So your expression is not quite correct. (April 1/08)
In 2008 assignment 5
question 5 (problem 10.2 of the book), what is the relationship between the
current source (is) and the load current (iL)?
I wanted to consider the whole thing as a serial circuit but then could not
link iL to is. Maybe there
is a lack of info in the question.
The electrical load torque Tg = 1/alpha_r x i_L, where alpha_r = Cr i_s (t). So the field current i_s affects the coupling coefficient between the torque into the generator and the current i_L in the electrical circuit. In the mixed model problems in the course text, not all the necessary information is stated explicitly in the problem. You may have to look up some of the constitutive relationships in other parts of the text. (Mar. 30/08)
Regarding Q5 - last
part: if I take a derivative of a differential equation that contains “i_s”, does “d/dt(i_s)” appear or we can assume it
as constant for an specific input value? (I mean that deviation from nominal
value remains constant through the actual operation or simulation. i.e. here by
(i_s) I mean (i_s)_hat)
There is nothing in the problem to tell us that the input i_s is constant so we need both i_s “bar” (normal operating point) and i_s “hat” (linear perturbation). We can assume, however, that second-order terms can be neglected. (April 1/08)
I’m not sure how to
do the Runge-Kutta integration. Do I need to have the
equations in first-order form?
Yes. In assignment 4, we do several problems in which we restate the governing equations in first-order form. To do this, we typically have to assign some new variables. This new set of equations has a common form: the derivative of one of our new variables is equal to a function of some (or all) of the new variables, as well as time. This arrangement allows us to use a solver such as Runge-Kutta to find an estimate of a variable at a discrete time in the future. The solver does this by extrapolating from the current time t by a time interval h to find an estimate for the variable at the time, using the expression for the derivative of that variable. Note that the expression for the derivative may be in terms of any or all of the variables; and so these variables will have to be found at the time of interest. This means that we can do a set of extrapolations at each incremental advance in time, thereby solving for all of the variables being integrated over time. If we have two first-order equations of motion, we'll do two separate solutions to find the values at t+h for each of the two variables. Then we do the next discrete time point at t+h+h, and so on, until we’ve produced a series of solutions for the set of variables for a series of discrete points in time within the time interval of interest. There is a description of how to do RK integration in the course text, and a worked example is posted in the Examples directory on the course web site. For a solver such as ode45, Matlab requires a function to evaluate the set of derivatives of the variables of interest. The handout for Lab 7 has an explanation of how to write such a function. (March 17/08)
I don’t know how to
start the Picard method integration in 2008 Assignment 3 question 2.
Recall that Picard’s method is an iterative approach to estimating the value of a function at a later time t, where we know the derivative and the initial values at time t=0. The system of interest is second-order, and so we use the first-order form, which has two equations. Each of these equations is generally in terms of all the variables of the system, say y1 and y2, as well as time. Picard’s method is a way to generate a power series with respect to time that gives an increasingly accurate expression of how the actual system behaves between time t=0 and t. So now, if we have two derivatives, say f1(y1,y2,t) and f2(y1,y2,t), we use the previous estimate of the expressions for y1 and y2 to improve the expression for each (adding more terms to the series which is a function of time). So, to reiterate, if we have two first-order equations, we should have two expressions to evaluate. We generate new expressions with respect for time for each of the two variables by respectively integrating the expressions for the derivatives (which are generally functions of one or maybe even both of the two variables, where the expression for each variable is taken from the previous step). This process is similar to doing iteration by total steps for equilibrium systems.
(Mar. 3/08)
I have a question regarding the course notes. In the notes for lecture 4, there
is an example regarding the buckling of a two-bar linkage joined with springs.
I wonder why you assume that the reaction forces are along the bars. It is done
before assuming small angle and even small angle assumption does not result in
such assumption.
This analysis approach for buckling assumes thin rigid bodies, with forces acting only axially through the linkages and at the nodes. This is a simplification of the physics of how forces flow in the links. You are right that the small angle assumption only simplifies the equations by limiting the geometry (the range of the angles) in the problem. (Mar 3/08) (Note for 2010: lecture 4 in 2010 is not the same: eigenvalue problems are introduced in Lecture 7.)
Why do we have to use
the general formulation for some of these problems? The mesh current analysis
method I learned in my undergrad electric circuits
course seems more intuitive.
This is a very good question. Our solution method is more generally applicable to a broad range of systems, and can handle nonlinearities in a particular system. In this way, we learn how to model systems in general, not particular methods for special kinds of systems. Many solution methods rely on superposition, and so they are only valid for linear systems. The systems we study are pretty simple, but the approach scales up to large - even nonlinear - systems. Note that on assignments and tests it is important to show your understanding of our general method, not just getting an answer by any means. (Feb 23/08)
I have a question
related to Eigenvalue Problems. In the electrical
problem, we use v = v'e^(iwt), i = i'
e^(iwt) to form the governing equation;
however, it seems we use u = u sinwt to
deal with other types of problems except electrical problems. Is that
true? Thank you for your time.
In these type of eigenvalue questions, a harmonic trial function is used that is differentiable and yields an equation in terms of the original trial function (when differentiated twice). In both cases, the resulting equations reveal the change of sign between element types that represents the physical situation of energy exchange between elements in a critical loading scenario.
For a linear eigenvalue problem, Xsin(wt) is a sufficiently useful trial function, because there is neither a forcing function to change the energy in the system (adding or removing) nor any energy dissipating elements in the system. Xsin(wt) is less generally useful than e^(iwt). Recall the Euler identity
e^(iwt) = cos(wt) + i sin(wt)
so e^(iwt)
is a more general form for sinusoidal inputs. When using e^(iwt) the original trial function appears in each
differentiation (not just reappearing in the second derivative). This is the
basis for the
What’s the algorithm
for Cholesky decomposition, and how do I implement it
in Matlab?
Cholesky decomposition is a variation of Gaussian elimination. Wikipedia has an entry on its formulation. In Matlab, if you have a matrix [B], use the command chol(B) to find the upper triangular matrix form, which is the transpose of [L]. (Feb 14/08)
Do eigenvalues only appear on the diagonal?
It depends on the nature of the system. If there is only one A-type energy storage element in a governing equation (such as a mass in an oscillating mechanical system, or an inductance in an electrical system), then there will be one eigenvalue in that equation. In our course, we consider cases where the matrix [B] is symmetric and positive definite. Diagonal forms of [B] fall into that category. (Feb. 14/08)
Are eigenvectors always loop variables like
voltages and displacements?
Eigenvectors can be expressed for node variables as well. For example, in a buckling problem, the forces may be of interest, in which case the eigenvectors would be in terms of node variables. (Feb. 11/08)
In 2008 Assignment 2
question 3, are the natural frequencies just the eigenvalues?
Also, how do you represent the mode shapes for an electrical system?
For your first question, the short answer is no; however, the natural frequency is related to the eigenvalue. The eigenvalue λ (lambda) is a relationship amongst certain system parameters that affect the critical condition of the system. For the mechanical systems that we have considered in the course, mass m, spring constant k, and natural frequency ω (omega) are the parameters embedded in the eigenvalue. Since the damping coefficient b is zero, this system becomes a linear eigenvalue problem.
Here’s the slightly longer answer. When you use the harmonic function to generate governing equations with respect to the chosen variable type (getting rid of the integrals & derivatives…) you end up with a solution that includes the frequency at which the system will oscillate if it is disturbed in some way. The system will have a set of these frequencies, one for each solution (eigenvector). The way that the system can behave is limited to combinations of the eigenvectors. (All the sounds that come out a musical instrument are combinations of a finite number of frequency responses.) If there is a dissipative element in the system (like a mechanical damper or an electrical resistor), then the system now has a different solution than if the damper were not there; and the system response is a damped oscillation, which looks like a sinusoidal wave that peters out with time.
The mode shapes for an electrical system will be harmonic responses, just like to the mechanical system with carts that we examined in the lecture. In the case of an electrical system, the eigenvectors are of course expressed as voltage loop variables rather than displacements. The solution for lambda = 3 has a shorter wavelength than that of lambda =1/2 because the natural frequencies are different.
Analogies between different types of systems used to be really important for doing simulations, before digital computers were fast enough to present results in real-time. Analog computers used electrical circuits that represented the behaviour of mechanical systems (or other types of systems), and plotted an output that corresponded to the output expected from the system being modeled. Which was pretty cool, way back in the 20th century. (Jan. 31/08)
I am having some difficulties with the first question of 2008 Assignment 2. Part (a) worked out fine but part (b) is diverging for some reason. Is that supposed to happen, or am I doing something wrong?
If you check the row-sum criterion, I think you’ll find that there is no guarantee of convergence in this case. So it is entirely possible that an iterative method might diverge. The row-sum criterion only tells you when convergence is guaranteed. If any of the alphas is equal to one or greater than one, the system may still converge; but it does not tell you whether a particular method will converge or diverge. If the estimates keep oscillating or begin to diverge after a couple of iterations, then the solution method isn't working and you should either change the initial guess or try another method. Although the Gauss-Siedel method has the advantage of faster convergence, it can also diverge quickly too. (Jan. 25/08, with update Feb 2/08)
How do we non-dimensionalise a variable?
When non-dimensionalising, it is important to remember that
we are trying to get an expression that has no units. The coefficient constant
on the right-hand side of the governing equation has an engineering unit
associated with it (volts,
Does the node variable change as it goes through an element?
No, it doesn’t. What goes into the element has to come out. For example, no current i is lost going through a resistor; what happens is that the loop variable acting across the element is affected (as a function of the constitutive relationship term R). So the voltage V across the resistor changes as i = (1/R) Vi. Note that the constitutive relationship is expressed with the loop variable on the x -axis; this is because energy is a function of loop variables.
Mechanical systems have force as the node variable, so the force that goes into the node (e.g. a spring) is the same as what will come out at the other end. The displacement u changes with force F as a function of the constitutive relationship k according to F = k u for a linear spring). (Jan. 10/08)
What’s the difference between an eigenvalue
problem and a propagation problem?
Both types of problem consider the dynamic behaviour of a system. In the case of a propagation problem, the system typically has non-conservative forces acting on the system, which means that energy can enter or leave the system over time. The state variables are solved as a function of time. (State variables are the variables in which the governing equations are formulated, either using loop variables or node variables.) In the case of an eigenvalue problem, the system conserves energy. The problem is formulated as a set of differential equations in homogeneous form, which means that there is no forcing function. If the case of a linear eigenvalue system, if there is a momentary disturbance to the system, then it will oscillate in a finite number of ways, moving energy amongst elements, but neither gaining nor losing energy in the overall system. (Jan. 10/08)