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).

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 (m_{3}) and wobbling (m_{4}) masses,
and the other two represent the lower body rigid (m_{1}) and wobbling (m_{2})
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 m_{1} is in direct contact with
the ground (so x_{1}(t)=0). The simplified model is
redrawn below.

a) Derive equations of motion.

b) Derive a state space model with three outputs, x_{2}, x_{3}
and x_{4}. Note there is no input, so you don't need a **
B** or **D** matrix.

c) The paper gives model values as m_{1}=6.15 kg, m_{2}=6 kg, m_{3}=12.6
kg, m_{4}=50.3 kg, k_{1}=k_{2}=6000 N/m, k_{3}=k_{4}=10000
N/m, k_{5}=18000 N/m, b_{1}=300 kg/s, b_{2}=650 kg/s
and b_{4}=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 v_{2}(0)=0.96 m/s, and v_{3}(0)=v_{4}(0)=2.0
m/s. Assume initial positions are all zero. Plot x_{2}(t),
x_{3}(t) and x_{4}(t). Matlab's "initial"
command may be used.

a) the two spring attached to mass 4 act as one larger spring so let k_{4}+k_{5}=k_{45}=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)

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:

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.

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").

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

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

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

Given the system:

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

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

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.

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, C_{1}=3000 J/°K.
The inner region (Region 2) has a thermal capacitance, C_{2}=1000 J/°K.
The thermal resistance between the external environment and Region 1 is R_{1a}=2°K/W.
The resistance between the inner and outer regions is R_{12}=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

- Draw an equivalent electrical circuit that represents:
- 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.
- 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).
- Develop a state space model (with two outputs - the temperature of each region) for
- the bolus in a conventional oven. The input is the external temperature.
- the bolus in a microwave oven. The input is the power from the microwave.
- If the bolus is initially at ambient temperature (call this 0°) and is suddenly placed in a conventional oven at 200 °K above ambient:
- plot the temperature of the two regions.
- approximately how long does it take region 2 to reach 50° above ambient?
- If the bolus is initially at ambient temperature (call this 0°) and is suddenly placed in a microwave oven that generate 1500 Watts:
- plot the temperature of the two regions.
- approximately how long does it take region 2 to reach 50° above ambient?

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.

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 (*z*ero *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) e_{s}(t)
represents the temperature pulse. R_{1} 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 (e_{s}(t)-e_{1}(t))/R_{1})).
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 R_{2} 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 R_{4}.

**Analogous statements about the electrical system:**

For the system in this problem (see image above) e_{s}(t)
represents ** a voltage** pulse. R

_{1}is the

**resistance (resistance to**

*electrical***flow) between the surface of the paper (where the pulse occurs) and the yellow crystals (the**

*current***flow is proportional to the**

*current***difference divided by the resistance (e**

*voltage*_{s}(t)-e

_{1}(t))/R

_{1})). As

**flows into the yellow crystals their temperature increases over time**

*current**(just as the voltage of a capacitor increases as current flows into it)*. As the yellow crystals'

**increases, some**

*voltage***will flow into the magenta crystals through**

*current***resistance R**

*an electrical*_{2}between the two sets of crystals. As the magenta crystals

**increases, some**

*voltage***will flow into the cyan crystals causing their**

*current***to rise. Finally some of the**

*voltage***current**from the cyan crystals can flow into the paper through R

_{4}.

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 C_{1} is e_{1}(t), the voltage across C_{2}
is e_{2}(t) and the voltage across C_{3} is e_{3}(t).

The input is e_{s}(t). Choose as state variables e_{1}(t),
e_{2}(t) and e_{3}(t). For now consider just one output, e_{1}(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)=(s**I**-**A**)^{-1},
and use it to find the transfer function E_{1}(s)/E_{s}(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, e_{1}(t), e_{2}(t)
and e_{3}(t) and plot the unit step response of the system (there should
be three outputs).

Note that, as expected, e_{1}(t)
changes most quickly and has the highest steady state value, and e_{3}(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/

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)

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 e_{s}(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 e_{1}(t) to reach 200 V, t_{i1}. Clearly
explain your method.

2) estimate how long it takes e_{2}(t) to reach 150 V, t_{i2}.

3) estimate how long it takes e_{3}(t) to reach 100 V, t_{i3}.

**ii) **If e_{s}(t) is a 190 V step input:

a)estimate how long it takes e_{2}(t) to reach 150 V, t_{ii2}.

b)estimate how long it takes e_{3}(t)to reach 100 V, t_{ii3}.

**iii) **If e_{s}(t) is a 140 V step, estimate how long it takes e_{3}(t)
to reach 100 V, t_{iii3}.

Note,

- t
_{i1}should be less than t_{i2}or t_{i3}. - t
_{ii2}should be less than t_{ii3}.

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

Just for fun (nothing to solve)

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

Answer