E12 Assignment 14

Description: 

For this assignment you may use MATLAB wherever you would like.

You solved for the FBD and SS model of problem 1 last week -  you can simplly use the results this week (answers are posted online).

Problems: 
1. Body with wobbling masses

The system shown below is from an article entitled "Effects of Muscle Fatigue on the Ground Reaction Force and Soft-Tissue Vibrations During Running: A Model Study," by Nikooyan and Zadpoor in IEEE Transactions on Biomedical Engineering, Vol  59, No. 3, page 797, March 2012.  The model is that of the human body as it comes in contact with the ground.  The paper examines fatigue as indicated by the amount of shaking of m2 and m4.


The paper describes the model as follows "Two of these masses represent the upper body rigid (m3) and wobbling (m4) masses, and the other two represent the lower body rigid (m1) and wobbling (m2) masses.  The springs and dampers represent the stiffness and damping properties of the human body hard and soft tissues."  As the foot hits the ground, there is an opposing force determined by the "ground reaction model" that depends in a complex way on the kind of turf and the kind of shoes being worn, and is too complex to consider here.

We will simplify the model by assuming that m1 is in direct contact with the ground (so x1(t)=0).   The simplified model is redrawn below.

a) Derive equations of motion.

b) Derive a state space model with three outputs, x2, x3 and x4.  Note there is no input, so you don't need a B or D matrix.

c) The paper gives model values as m1=6.15 kg, m2=6 kg, m3=12.6 kg, m4=50.3 kg, k1=k2=6000 N/m, k3=k4=10000 N/m, k5=18000 N/m, b1=300 kg/s, b2=650 kg/s and b4=1900 kg/s.   Examine the roots of the characteristic equation (Matlab's "damp", "zpkdata" or "pzmap" can be used) and state whether or not you think the system will exhibit oscillatory behavior.

d) During a run, just as the foot strikes the ground, initial conditions are given in the paper for the velocities v2(0)=0.96 m/s, and v3(0)=v4(0)=2.0 m/s.  Assume initial positions are all zero.  Plot x2(t), x3(t) and x4(t).  Matlab's "initial" command may be used.

Solution: 

a) the two spring attached to mass 4 act as one larger spring so let k4+k5=k45=28000 N/m.

b)

c)  Some of the roots of the characteristic equation are complex, so we expect some oscillation (ζ=0.33).

%% Matlab commands
% Define constants
m2=6; m3=12.6; m4=50.3;
k1=6000; k2=6000; k3=10000; k45=28000;
b1=300; b2=650; b4=1900;

% Define state matrices
A = [[0 1 0  0 0 0];
    [-k2-k3 -b2 k3  0 0 0]/m2;
    [0 0 0  1 0 0];
    [k3 0 -k1-k3-k45 -b1-b4 k45 b4]/m3;
    [0 0 0  0 0 1];
    [0 0 k45  b4 -k45 -b4]/m4];
B=[0 0 0 0 0 0]';
C=[1 0 0 0 0 0;
    0 0 1 0 0 0;
    0 0 0 0 1 0];
D=0;

sys=ss(A,B,C,D);
damp(sys)
 %% Matlab output
   Eigenvalue            Damping     Freq. (rad/s)  
                                                         
 -3.91e+000 + 1.11e+001i     3.32e-001      1.18e+001    
 -3.91e+000 - 1.11e+001i     3.32e-001      1.18e+001    
 -1.52e+001                  1.00e+000      1.52e+001    
 -4.39e+001                  1.00e+000      4.39e+001    
 -6.57e+001                  1.00e+000      6.57e+001    
 -1.88e+002                  1.00e+000      1.88e+002  

d)

q0=[0 0.96 0 2.0 0 2.0]';
initial(sys,q0)

2. State Space to Xfer Func - 2 systems

a) Find the transfer function of the system defined by:

Recall that the inverse of a 2x2 matrix is given by:

b) Find the transfer function of the system defined by:

Solution: 

a)

>> syms s
>> A=[-5 -2; 2 0];  B=[2 1]';  C=[3 2];  D=0;

>> s*eye(2)-A
ans =
[ s + 5, 2]
[    -2, s]

>> Phi=inv(ans)
Phi =
[ s/(s^2 + 5*s + 4),      -2/(s^2 + 5*s + 4)]
[ 2/(s^2 + 5*s + 4), (s + 5)/(s^2 + 5*s + 4)]

>> H=simple(C*Phi*B+D)
H =
(8*s + 12)/(s^2 + 5*s + 4)

b)

>> A=[0 2; -2 -5];  B=[1 2]';  C=[2 3];  D=0;
>> H=simple(C*inv(s*eye(2)-A)*B+D)
H =
(8*s + 12)/(s^2 + 5*s + 4)

Note: Systems are the same, but state variables are switched.

3. Xfer Function to State Space, 2nd order, CCF, OCF

For the transfer function

a) Find the state space representation in Controllable Canonic form. CCF form is described here). You may simply use the stated result - there is no need to rederive anything.

b) Define the state space system in MATLAB (with the "ss" command) and then use it to verify your answer by finding the transfer function ("tfdata").

c) Repeat for Observable Canonic form. OCF form is described here). You may simply use the stated result - there is no need to rederive anything.

d) Check with MATLAB (as in part "b").

Solution: 

a)    

 b)        Check with Matlab

>> sys=ss([0 1;-6 -5],[0 1]', [1 1], 0);
>> [n,d]=tfdata(sys,'v')
n =         0    1.0000    1.0000
d =    1.0000    5.0000    6.0000

c) 

d)

>> sys=ss([-5 1;-6 0],[1 1]', [1 0], 0); 
>> [n,d]=tfdata(sys,'v') 
n = 0      1.0000 1.0000 
d = 1      5      6 
4. Xfer function to State Space, 3rd order, CCF

Find the controllable canonic form of a system described by the transfer function (CCF form is described here)

Solution: 

Check with MATLAB

>> sys=ss([0 1 0; 0 0 1; -1 -1 -2],[0 0 1]',[-2 -2 -6],3);

>> [n,d]=tfdata(sys,'v')
n =
    3.0000   -0.0000    1.0000    1.0000
d =
    1.0000    2.0000    1.0000    1.0000
5. State Space to State Space, 2nd order

Given the system:

a) Transform (by hand) to a new set of state variables where,

b) Check with MATLAB using the "ss2ss" command.

Solution: 

a)

b)

>> sys=ss([0 2;-2 -5],[1 2]',[2 3],0)
sys =
  a = 
       x1  x2
   x1   0   2
   x2  -2  -5
  b = 
       u1
   x1   1
   x2   2
  c = 
       x1  x2
   y1   2   3
  d = 
       u1
   y1   0
Continuous-time state-space model.

>> sysnew=ss2ss(sys,[0 1; 1 0])
sysnew =
  a = 
       x1  x2
   x1  -5  -2
   x2   2   0
  b = 
       u1
   x1   2
   x2   1
  c = 
       x1  x2
   y1   3   2
  d = 
       u1
   y1   0
 
Continuous-time state-space model.

Note: Systems are the same, but state variables are switched.

6. Heating Food in a Microwave

The next problems investigate why microwave ovens are faster than conventional ovens using a simplistic, but not altogether unreasonable, model.  The task is to heat a bolus of meat about the size and shape of an orange; we will model the bolus as two masses as shown.  The outer region (Region 1) has a thermal capacitance, C1=3000 J/°K.  The inner region (Region 2) has a thermal capacitance, C2=1000 J/°K.  The thermal resistance between the external environment and Region 1 is R1a=2°K/W.  The resistance between the inner and outer regions is R12=0.1°K/W.   A diagram is shown.

Note: the mechanism of heat transfer in a microwave oven is radiation - but we can assume this is just a constant power input

  1. Draw an equivalent electrical circuit that represents:
    1. the bolus in a conventional oven in which the temperature outside the bolus is controlled.  We will assume the temperature outside the bolus starts at ambient temperature (we consider this 0°).  When we put the food in the oven, we will model this as a step change in external temperature.
    2. the bolus in a microwave oven in which power is deposited directly in region 1 (but not in region 2 because the depth of penetration of microwaves in tissue is limited to 1-2 cm).  Assume the temperature in the oven is constant and equal to ambient temperature.  When we turn the microwave on it immediately begins depositing power in region 1 (but the air inside the microwave does not heat up appreciably).
  2. Develop a state space model (with two outputs - the temperature of each region) for
    1. the bolus in a conventional oven.  The input is the external temperature.
    2. the bolus in a microwave oven.  The input is the power from the microwave.
  3. If the bolus is initially at ambient temperature (call this 0°) and is suddenly placed in a conventional oven at 200 °K above ambient:
    1. plot the temperature of the two regions.
    2. approximately how long does it take region 2 to reach 50° above ambient?
  4. If the bolus is initially at ambient temperature (call this 0°) and is suddenly placed in a microwave oven that generate 1500 Watts:
    1. plot the temperature of the two regions.
    2. approximately how long does it take region 2 to reach 50° above ambient?
Solution: 

a)  Diagram on left shows conventional oven (oven temperature is input), diagram at the right shows a microwave oven (input power is input).

b) Note that the state space systems are the same except for the B matrix.  This is because the systems are identical with zero input (the voltage source becomes a short circuit, and the current source becomes an open circuit).

Convential Oven   Microwave Oven
Energy Balance
a=0)
 
  State Space

c)
>> C1=3000; C2=1000; R12=0.1;  R1a=2;
>> A=[-(1/R1a+1/R12)/C1 1/R12/C1; 1/R12/C2 -1/R12/C2];
>> B=[1/R1a/C1 0]';  C=[1 0; 0 1;];  D=0;
>> sys=ss(A,B,C,D);
>> sys.OutputName={'Temp 1','Temp 2'}; sys.InputName='Input TempPower';
>> sys.StateName={'T1','T2'};
>> step(200*sys)    % Note step size is 200.

 

It takes about 2380 seconds (≈40 minutes) to reach 50 degrees.

d) We need only change the "B" matrix and the amplitude of the input

>> sys.b=[1/C1 0]'; 
>> step(1500*sys)

Zoom in for detail.  It takes about 205 seconds (≈3.5 minutes) to reach 50 degrees.

7. Zink Printer

You can use MATLAB wherever you would like, but include any code that you use, and make sure it is well commented.


(you can skip this entire section if desired and move on to the problem)

Background

In the fall of 2009 a company (zink) debuted a color printing technology that requires no ink cartridges (zero ink = zink) to be used with portable devices like cameras.  As part of the development process they developed three chemicals that were normally colorless, but would change color when they reached a certain temperature.  One chemical would turn yellow at about 200°C, another would turn magenta at 150°C, and the third would turn cyan at 100°C.

The inks are then layered with yellow on top, magenta in the middle and cyan at the bottom. 

  • A very short intense pulse of heat would heat up the top layer enough to turn it yellow, but the pulse is short enough that the other two layers don't change color. 
  • A lower temperature (below 200°C) pulse would change the middle layer but not the top layer (because it was below 200°C) and not the bottom layer (because the pulse is short enough that the temperature remains below 100°C.
  • An even lower temperature pulse only effects the cyan layer.


image from IEEE web page, reference at end of this web page

Thermal systems can be modeled as electrical systems with voltage analogous to temperature and heat flow analogous to current. 

The model

Statements about the thermal system:

For the system in this problem (see image above) es(t) represents the temperature pulse.  R1 is the thermal resistance (resistance to heat flow) between the surface of the paper (where the pulse occurs) and the yellow (the heat flow is proportional to the temperature difference divided by the resistance (es(t)-e1(t))/R1)).  As heat flows into the yellow crystals their temperature increases over time.  As the yellow crystals' temperature increases, some heat will flow into the magenta crystals through a thermal resistance R2 between the two sets of crystals.  As the magenta crystals temperature increases, some heat will flow into the cyan crystals causing their temperature to rise.  Finally some of the heat from the cyan crystals can flow into the paper through R4.

Analogous statements about the electrical system:

For the system in this problem (see image above) es(t) represents a voltage pulse.  R1 is the electrical resistance (resistance to current flow) between the surface of the paper (where the pulse occurs) and the yellow crystals (the current flow is proportional to the voltage difference divided by the resistance (es(t)-e1(t))/R1)).  As current flows into the yellow crystals their temperature increases over time (just as the voltage of a capacitor increases as current flows into it).  As the yellow crystals' voltage increases, some current will flow into the magenta crystals through an electrical resistance R2 between the two sets of crystals.  As the magenta crystals voltage increases, some current will flow into the cyan crystals causing their voltage to rise.  Finally some of the current from the cyan crystals can flow into the paper through R4.

Caveat: the time constants associated with the thermal system are very short (but not made publically available by zink).  I have made the time constants of the problem much larger to make the numbers easier to work with.


The Problem

Consider the circuit shown all voltages are measured relative to ground.  The voltage across C1 is e1(t), the voltage across C2 is e2(t) and the voltage across C3 is e3(t).

The input is es(t).  Choose as state variables e1(t), e2(t) and e3(t).  For now consider just one output, e1(t).

a) Derive a state space representation in terms of the resistors and capacitors in the circuit.  What are A, b, C and d?

For a certain combination of resistors and capacitors the state space representation becomes:

 

b) Using these values, find an expression for Φ(s)=(sI-A)-1, and use it to find the transfer function E1(s)/Es(s). (Use a symbolic solver - you needn't do this by hand).

c) Verify your result from part b using MATLAB to find the transfer function using either the "ss2tf" function or the "tfdata" function. 

d) This is a third order system with real roots, find the three time constants associated with the system (use "damp" or "poles" or simply find the roots of the denominator of the transfer function.

e) Modify your state space representation so there are three outputs, e1(t), e2(t) and e3(t) and plot the unit step response of the system (there should be three outputs). 
Note that, as expected, e1(t) changes most quickly and has the highest steady state value, and e3(t) changes most slowly and has the lowest steady state value.

If you would like to see more information about zink, see: http://spectrum.ieee.org/consumer-electronics/gadgets/zink-inkless-printing-with-colorless-color/0 or http://www.zinkless.com/

Solution: 

a) Find each capacitor current (C·de/dt) as a sum of currents through resistors, and divide by C to get first order differential equations.

b)

R1=1; R2=2; R3=6; R4=60;
C1=0.1/12; C2=0.2/12; C3=0.4/12;

A=[ [-1/R1-1/R2 1/R2 0]/C1;,...
    [1/R2 -1/R2-1/R3 1/R3]/C2;...
    [0 1/R3 -1/R3-1/R4]/C3];
B=[1/R1/C1 0 0]';
C=[1 0 0];
D=0;
mySys=ss(A,B,C,D)

syms s
PHI=inv(s*eye(3)-A);
disp("Xfer func");
pretty(simple(PHI))
mySys =
  a = 
         x1    x2    x3
   x1  -180    60     0
   x2    30   -40    10
   x3     0     5  -5.5
  b = 
        u1
   x1  120
   x2    0
   x3    0
  c = 
       x1  x2  x3
   y1   1   0   0
  d = 
       u1
   y1   0
Continuous-time state-space model.


/    2                                                           \
| 2 s  + 91 s + 340      (2 s + 11) 60              1200         |
| -----------------,     -------------,             ----         |
|         #1                   #1                    #1          |
|                                                                |
|   (2 s + 11) 30    (2 s + 11) (s + 180)       (s + 180) 20     |
|   -------------,   --------------------,      ------------     |
|         #1                  #1                     #1          |
|                                                                |
|                                            2                   |
|        300             (s + 180) 10      (s  + 220 s + 5400) 2 |
|        ---,            ------------,     --------------------- |
\         #1                  #1                     #1          /

where
            3        2
   #1 == 2 s  + 451 s  + 13120 s + 41400
tfunc=simple(C*PHI*B+D);
pretty(tfunc)

        2
    (2 s  + 91 s + 340) 120
-------------------------------
   3        2
2 s  + 451 s  + 13120 s + 41400

c)

[n,d]=tfdata(mySys,'v')

% multiply n and d by 2 to compare with previous (since expression above has 2 with s^3).
n*2
d*2

n =
   1.0e+04 *
         0    0.0120    0.5460    2.0400
d =
   1.0e+04 *
    0.0001    0.0226    0.6560    2.0700
ans =
   1.0e+04 *
         0    0.0240    1.0920    4.0800
ans =
   1.0e+04 *
    0.0002    0.0451    1.3120    4.1400

d)

disp('Time constants:');
disp(-1./roots(d));
Time constants:
0.0052
0.0333
0.2784

e)

mySys.C=[1 0 0; 0 1 0; 0 0 1];
step(mySys)

8. Zink printer - quantitative results

Answer the following. You may use the results of the previous questions (read values off of the graph as accurately as you can (you may need to zoom in), or write code to find the answer- there is no need to do extensive calculations). You can use MATLAB wherever you would like, but include any code that you use, and make sure it is well commented.

i) If es(t) is a 250 V step input (Note - this problem is stated in terms of volts, but voltage is simply an analog for temperature):
     1) estimate how long it takes e1(t) to reach 200 V, ti1.  Clearly explain your method.
     2) estimate how long it takes e2(t) to reach 150 V, ti2.
     3) estimate how long it takes e3(t) to reach 100 V, ti3.

ii) If es(t) is a 190 V step input:  
     a)estimate how long it takes e2(t) to reach 150 V, tii2.
     b)estimate how long it takes e3(t)to reach 100 V, tii3.

iii) If es(t) is a 140 V step, estimate how long it takes e3(t) to reach 100 V, tiii3.

Note,

  • ti1 should be less than ti2 or ti3.
  • tii2 should be less than tii3.
Solution: 

Note: we can use the unit step response results, and scale to input voltage. For example, if step is 250 V then 200 V is scaled to 200/250=0.8.

%% E12 Zink2

R1=1; R2=2; R3=6; R4=60;
C1=0.1/12; C2=0.2/12; C3=0.4/12;

A=[ [-1/R1-1/R2 1/R2 0]/C1;,...
    [1/R2 -1/R2-1/R3 1/R3]/C2;...
    [0 1/R3 -1/R3-1/R4]/C3];
B=[1/R1/C1 0 0]';
C=[1 0 0; 0 1 0; 0 0 1];
D=0;
mySys=ss(A,B,C,D);

t=linspace(0,1,1001); % Sample output every msecond
[y]=step(mySys,t);

% Find first time value where temperature is above 200 degrees...
fprintf('\n\n\n');
fprintf('Part i,1   = %5.4f\n',min(t(y(:,1)>(200/250))));
fprintf('Part i,2   = %5.4f\n',min(t(y(:,2)>(150/250))));
fprintf('Part i,3   = %5.4f\n\n',min(t(y(:,3)>(100/250))));
fprintf('Part ii,1  = %5.4f\n',min(t(y(:,2)>(150/190))));
fprintf('Part ii,2  = %5.4f\n\n',min(t(y(:,3)>(100/190))));
fprintf('Part iii,1 = %5.4f\n',min(t(y(:,3)>(100/140))));
Part i,1   = 0.0440
Part i,2   = 0.0780
Part i,3   = 0.2130

Part ii,1  = 0.2330
Part ii,2  = 0.3000

Part iii,1 = 0.5210

% Or, even better, interpolate
fprintf('\n\n\n');
fprintf('Part i,1   = %5.4f\n',interp1(y(:,1),t,200/250));
fprintf('Part i,2   = %5.4f\n',interp1(y(:,2),t,150/250));
fprintf('Part i,3   = %5.4f\n\n',interp1(y(:,3),t,100/250));
fprintf('Part ii,1  = %5.4f\n',interp1(y(:,2),t,150/190));
fprintf('Part ii,2  = %5.4f\n\n',interp1(y(:,3),t,100/190));
fprintf('Part iii,1 = %5.4f\n',interp1(y(:,3),t,100/140));
Part i,1   = 0.0438
Part i,2   = 0.0779
Part i,3   = 0.2121

Part ii,1  = 0.2329
Part ii,2  = 0.2995

Part iii,1 = 0.5203
9. Hair Dryer in a box

Just for fun  (nothing to solve)

I used numbers for heat transfer between a steel box and air and got roughly the same results.

Solution: 

Answer