System Dynamics
System Dynamics
Home

Problems

Problem 14.1 (SIGMUND)

Consider a nonlinear capacitor with constitutive equation relating charge \(q_C\) and voltage \(v_C\): \[\begin{aligned} q_C = k v_C^{3/2} \end{aligned}\] with \(k\) a positive constant.

  1. Derive an elemental equation relating \(\diff v_C/ \diff t\) and \(i_C\) for the nonlinear capacitor.

  2. From the elemental equation, what is the voltage-dependent capacitance \(C(v_C)\)?

  3. Consider the RLC-circuit of figure 14.7, which includes the nonlinear capacitor. Derive a nonlinear state-space equation with state vector \[\begin{aligned} \bm{x} &= \begin{bmatrix} v_C & i_L \end{bmatrix}^\top. \end{aligned}\]

  4. For a constant input \(V_S(t) = 5\,\)V, derive the equilibrium state.

  5. Linearize the state-space equation about the operating point \[\begin{aligned} \bm{x}_o, \bm{u}_o &= \begin{bmatrix} 5\,\text{V} & 0\,\text{A} \end{bmatrix}^\top, \begin{bmatrix} 5\,\text{V} \end{bmatrix}. \end{aligned}\] Define the state equation matrices \(A\) and \(B\), the linearized state and input vectors \(\bm{x}^*\) and \(\bm{u}^*\), and the linearized state equation.

Circuit for problem 14.1.
Figure 14.7: Circuit for problem 14.1.

  1. The elemental equation relates \(i_C\) and \(v_C\), so time-differentiating the constitutive equation given yields \[\begin{aligned} i_C &= \frac{3 k} {2} v_C^{1/2}\frac{\diff v_C} {\diff t} \Rightarrow \\ \frac{\diff v_C} {\diff t} &= \frac{2} {3 k \sqrt{v_C}} i_C, \end{aligned}\] which is the (nonlinear) elemental equation.

  2. The linear elemental equation with capacitance \(C\) is \[\begin{aligned} \frac{\diff v_C} {\diff t} &= \frac{1} {C} i_C. \end{aligned}\] Comparing this to the nonlinear elemental equation from (a), the voltage-dependent capacitance \(C(v_C)\) of the nonlinear equation is related as \[\begin{aligned} \frac{1} {C(v_C)} &= \frac{2} {3 k \sqrt{v_C}} \Rightarrow \\ C(v_C) &= \frac{3 k \sqrt{v_C}}{2}. \end{aligned}\]

  3. The linear graph is straightforwardly bijective to the circuit diagram. The normal tree includes the source, capacitor, and resistor. The state variables, then, are \(v_C\) and \(i_L\). The state and input vectors are \[\begin{aligned} \bm{x} = \begin{bmatrix} v_C & i_L \end{bmatrix}^\top \quad \text{and} \quad \bm{u} = \begin{bmatrix} V_S \end{bmatrix}. \end{aligned}\] The elemental, continuity, and compatibility equations are as follows.

    Substituting the constraint equations into the elemental equations yields the following.

    Substituting the \(R\) equation in the \(L\) equation to eliminate \(v_R\), the state equation emerges: \[\begin{aligned} \frac{\diff \bm{x}}{\diff t} &= \bm{f}(\bm{x},\bm{u}) \\ &= \begin{bmatrix} \frac{2} {3 k \sqrt{v_C}} i_L \\ \frac{1} {L} (-R i_L + V_S - v_C) \end{bmatrix}. \end{aligned}\]

  4. An equilibrium state occurs when \(\diff\bm{x}/\diff t = \bm{0}\). Substituting this and \(V_S = 5\) into the state equation, we obtain the nonlinear algebraic system \[\begin{aligned} \begin{bmatrix} 0 \\ 0 \end{bmatrix} = \begin{bmatrix} \frac{2} {3 k \sqrt{v_C}} i_L \\ \frac{1} {L} (-R i_L + 5 - v_C) \end{bmatrix} \end{aligned}\] If \(v_C \ne 0\), the first equation implies \(i_L = 0\) and the second equation further implies \(v_C = 5\) V. If \(v_C = 0\), the first equation implies \(i_L = 0\) but the second equation then states \(0 = 5\), a contradiction implying this is not a solution. Therefore, the equilibrium is \[\begin{aligned} \overline{\bm{x}}, \overline{\bm{u}} = \begin{bmatrix} 5 & 0 \end{bmatrix}^\top, \begin{bmatrix} 5 \end{bmatrix}. \end{aligned}\]

  5. The equilibrium state and input from above are the operating point about which we are linearizing the state equation. The linearized system state vector is \[\begin{aligned} \bm{x}^* &= \bm{x} - \bm{x}_o \tag{by definition}\\ &= \begin{bmatrix} v_C - 5 \\ i_L \end{bmatrix}. \end{aligned}\] The linearized system input vector is \[\begin{aligned} \bm{u}^* &= \bm{u} - \bm{u}_o \tag{by definition}\\ &= \begin{bmatrix} V_S - 5 \end{bmatrix}. \end{aligned}\] The Jacobian matrix for \(A\) is \[\begin{aligned} A &= \left.\begin{bmatrix} \frac{-1} {3 k v_C^{3/2}} i_L & \frac{2} {3 k \sqrt{v_C}} \\ -1/L & -R/L \end{bmatrix}\right|_{\bm{x}_o,\bm{u}_o} \\ &= \begin{bmatrix} 0 & \frac{2} {3 k \sqrt{5}} \\ -1/L & -R/L \end{bmatrix} \end{aligned}\] and the Jacobian for the \(B\) matrix is \[\begin{aligned} B &= \left.\begin{bmatrix} 0 \\ 1/L \end{bmatrix}\right|_{\bm{x}_o,\bm{u}_o} \\ &= \begin{bmatrix} 0 \\ 1/L \end{bmatrix}. \end{aligned}\]

Problem 14.2 (FREUD)

In problem 14.1, you derived a nonlinear state-space model for the RLC circuit of figure 14.7, which includes a nonlinear capacitor, and linearized the state equation about an operating point. Use these results to perform the following analysis.

  1. Write a program to simulate the nonlinear state-space model for initial condition \(\bm{x}(0)=\begin{bmatrix}1 & 0 \end{bmatrix}^\top\) and step input \(\bm{u}(t) = 5 u_s(t)\). Let \(R = 10~\Omega\), \(L = 1\) mH, and \(k = 10^{-6}\). Try simulating for \(1\) ms.

  2. Add to the program the simulation of the linearized system for the same initial condition and input.

  3. Compare (by graphing) the nonlinear and linearized step responses. (Don’t forget that \(\bm{x}^* \ne \bm{x}\)!)

clear; close
save_figures = 1;

Part a: nonlinear solution

We begin by defining the nonlinear state equation in the following function file (could be a local function).

type ode_15_1a.m
function f = ode_nonlin(t,x)
    k = 1e-6;
    R = 10;
    L = 1e-3;
    VS = 5;
    f = [ ...
      2/(3*k*sqrt(x(1)))*x(2); ...
      1/L*(-R*x(2) + VS - x(1)) ...
    ];
end

Now let’s use ode45 for the simulation.

x0 = [1;0]; % initial condition
t_a = linspace(0,1e-3,200);
[~,x_a] = ode45( ...
    @ode_15_1a, ... % ODE function handle
    t_a, ... % time array
    x0 ... % initial state
);

Now we plot. Two vertical axes are used so the full range of each state variable is usefully displayed.

figure
yyaxis left
plot(1e3*t_a,x_a(:,1),'linewidth',1.5)
xlabel('time (ms)')
ylabel('v_C (V)')

yyaxis right
plot(1e3*t_a,1e3*x_a(:,2),'linewidth',1.5)
ylabel('i_L (mA)')
Nonlinear step response.
Figure : Nonlinear step response.

Part b: linear solution

Instead of using the usual linear simulatino methods, for parallelism with the nonlinear solution, we begin by defining the linear state equation in the following function file (could be a local function).

type ode_15_1b.m
function f = ode_nonlin(t,x)
    % not performant really
    k = 1e-6;
    R = 10;
    L = 1e-3;
    VS = 5;
    u = [VS];
    A = [ ... 
      0, 2/(3*k*sqrt(5)); ...
      -1/L, -R/L ...
    ];
    B = [ ...
      0; 1/L ...
    ];
    f = A*x + B*u; % linear state eq!
end

Now let’s use ode45 for the simulation.

x0 = [1;0]; % initial condition
t_a = linspace(0,1e-3,200);
[~,x_b] = ode45( ...
    @ode_15_1b, ... % ODE function handle
    t_a, ... % time array
    x0 ... % initial state
);

We could plot the linear solution alone, but let’s just plot them together, which is the next part.

Part c: plotting the linear and nonlinear solutions together

We use a similar approach.

figure
yyaxis left
hold on
plot(1e3*t_a,x_a(:,1),'linewidth',1.5)
plot(1e3*t_a,x_b(:,1),'linewidth',1.5)
xlabel('time (ms)')
ylabel('v_C (V)')

yyaxis right
plot(1e3*t_a,1e3*x_a(:,2),'linewidth',1.5)
plot(1e3*t_a,1e3*x_b(:,2),'linewidth',1.5)
hold off
ylabel('i_L (mA)')
legend('nonlinear v_C', ...
    'linear v_C', ...
    'nonlinear i_L', ...
    'linear i_L' ...
)

These plots show that the linear solution is pretty good, but would not be sufficient for precise analysis.

Comparison of nonlinear and linear step responses.
Figure : Comparison of nonlinear and linear step responses.

Problem 14.3 (FRANZ)

A nonlinear diode model gives a diode’s elemental equation to be \[\begin{aligned} i_D = I_s (\exp{(v_D/V_\text{TH})} - 1). \end{aligned}\] We let the saturation current be \(I_s = 10^{-12}\) A and the thermal voltage be \(V_\text{TH} = 0.025\) V. Considering this nonlinear diode model for the circuit of figure 14.8.

Circuit for .
Figure 14.8: Circuit for .

  1. Derive a nonlinear state-space equation with state vector \[\begin{aligned} \bm{x} &= \begin{bmatrix} v_C & i_L \end{bmatrix}^\top. \end{aligned}\] Hint: include the diode in your normal tree.

  2. For a constant input \(V_S(t) = 0\,\)V, derive the equilibrium state.

  3. Linearize the state-space equation about the operating point \[\begin{aligned} \bm{x}_o, \bm{u}_o &= \begin{bmatrix} 0\,\text{V} & 0\,\text{A} \end{bmatrix}^\top, \begin{bmatrix} 0\,\text{V} \end{bmatrix}. \end{aligned}\] Hint: \(\diff{} \ln(x)/\diff{x} = 1/x\). Define the state equation matrices \(A\) and \(B\), the linearized state and input vectors \(\bm{x}^*\) and \(\bm{u}^*\), and the linearized state equation.

  1. The linear graph is straightforwardly bijective to the circuit diagram. The normal tree includes the source, capacitor, resistor, and diode. The state variables, then, are \(v_C\) and \(i_L\). The state and input vectors are \[\begin{aligned} \bm{x} = \begin{bmatrix} v_C & i_L \end{bmatrix}^\top \quad \text{and} \quad \bm{u} = \begin{bmatrix} V_S \end{bmatrix}. \end{aligned}\] The elemental, continuity, and compatibility equations are as follows.

    Substituting the constraint equations into the elemental equations yields the following.

    Solving the \(D\) equation for \(v_D\) gives \[\begin{aligned} v_D = V_\text{TH} \ln(i_L/I_s + 1). \end{aligned}\] Substituting this and the \(R\) equation into the \(L\) equation to eliminate \(v_D\) and \(v_R\), the state equation emerges: \[\begin{aligned} \frac{\diff \bm{x}}{\diff t} &= \bm{f}(\bm{x},\bm{u}) \\ &= \begin{bmatrix} \frac{1} {C} i_L \\ \frac{1} {L} (-V_\text{TH} \ln(i_L/I_s + 1) - R i_L + V_S - v_C) \end{bmatrix}. \end{aligned}\]

  2. An equilibrium state occurs when \(\diff\bm{x}/\diff t = \bm{0}\). Substituting this and \(V_S = 0\) into the state equation, we obtain the nonlinear algebraic system \[\begin{aligned} \begin{bmatrix} 0 \\ 0 \end{bmatrix} = \begin{bmatrix} \frac{1} {C} i_L \\ \frac{1} {L} (-V_\text{TH} \ln(i_L/I_s + 1) - R i_L + 0 - v_C) \end{bmatrix} \end{aligned}\] The first equation implies that for equilibrium \(i_L = 0\). The second equation, with \(i_L = 0\), implies \(v_C = 0\). Therefore, the equilibrium is \[\begin{aligned} \overline{\bm{x}}, \overline{\bm{u}} = \begin{bmatrix} 0 & 0 \end{bmatrix}^\top, \begin{bmatrix} 0 \end{bmatrix}. \end{aligned}\]

  3. The equilibrium state and input from above are the operating point about which we are linearizing the state equation. The linearized system state vector is \[\begin{aligned} \bm{x}^* &= \bm{x} - \bm{x}_o \tag{by definition}\\ &= \begin{bmatrix} v_C - 0 \\ i_L - 0 \end{bmatrix}\\ &= \begin{bmatrix} v_C \\ i_L \end{bmatrix}. \end{aligned}\] The linearized system input vector is \[\begin{aligned} \bm{u}^* &= \bm{u} - \bm{u}_o \tag{by definition}\\ &= \begin{bmatrix} V_S - 0 \end{bmatrix}\\ &= \begin{bmatrix} V_S \end{bmatrix}. \end{aligned}\] The Jacobian matrix for \(A\) is \[\begin{aligned} A &= \left.\begin{bmatrix} 0 & 1/C \\ -1/L & \frac{-V_\text{TH}/L}{i_L + I_s} - R/L \end{bmatrix}\right|_{\bm{x}_o,\bm{u}_o} \\ &= \begin{bmatrix} 0 & 1/C \\ -1/L & -(V_\text{TH}/I_s + R)/L \end{bmatrix} \end{aligned}\] and the Jacobian for the \(B\) matrix is \[\begin{aligned} B &= \left.\begin{bmatrix} 0 \\ 1/L \end{bmatrix}\right|_{\bm{x}_o,\bm{u}_o} \\ &= \begin{bmatrix} 0 \\ 1/L \end{bmatrix}. \end{aligned}\]

Problem 14.4 (KAFKA)

Let the nonlinear state equation of a circuit like figure 14.8, including a diode, be \[\begin{aligned} \frac{\diff \bm{x}}{\diff t} &= \bm{f}(\bm{x},\bm{u}) \\ &= \begin{bmatrix} \frac{1} {C} i_L \\ \frac{1} {L} (-V_\text{TH} \ln(i_L/I_s + 1) - R i_L + V_S - v_C) \end{bmatrix}. \end{aligned}\]

  1. Write a program to simulate the nonlinear state-space model for initial condition \(\bm{x}(0)=\begin{bmatrix}0 & 0 \end{bmatrix}^\top\) and input \(\bm{u}(t) = 1 + 0.1 \cos(8000 \pi t)\). Let \(I_s = 10^{-12}\) A, \(V_\text{TH} = 25\) mV, \(R = 10~\Omega\), \(L = 1\) mH, and \(C = 10 \mu\)F. Try simulating for \(1\) ms. Hint: the ode is stiff, so simulate with .

  2. Add to the program the simulation of the linearized system (with operating point \(\bm{x}_o = \begin{bmatrix}0 & I_s\end{bmatrix}^\top\), \(u_o = 0\)) with \(A\) and \(B\) matrices \[\begin{aligned} A = \begin{bmatrix} 0 & 1/C \\ -1/L & -(V_\text{TH}/(2 I_s) + R)/L \end{bmatrix} \quad \text{and} \quad B = \begin{bmatrix} 0 \\ 1/L \end{bmatrix} \end{aligned}\] for the same initial condition and input.

  3. Compare (by graphing) the nonlinear and linearized step responses. (Don’t forget that \(\bm{x}^* \ne \bm{x}\)!)

Begin by defining parameters globally.

global R L C VTH IS
R = 10; % Ohms
L = 1e-3; % H
C = 10e-6; % F
VTH = 25e-3; % V
IS = 1e-12; % A

Now define the initial condition, input, and a time array for simulation.

x0 = [0;0]; % initial state
u = @(t) 1 + .1*cos(4*2*pi*1e3*t);
t_a = linspace(0,1e-3,700);

Part a: nonlinear solution

We begin by defining the nonlinear state equation in the following function file (could be a local function).

type kafka_ode.m
function dxdt = kafka_ode(t,x,u)
    global R L C VTH IS
    dxdt = [ 1/C*x(2);
        1/L*(u(1) - R*x(2) - VTH*log(x(2)/IS+1) - x(1)) ];
end

Let’s first try ode45 for the simulation.

[~,x] = ode45( ...
    @(t,x) kafka_ode(t,x,u(t)), ... % ODE ... 
    t_a, ... % time array or span
    x0 ... % initial state
);

We now have our initial solution array x, which represents our state solution \(x(t)\). Let’s plot the solution.

figure
yyaxis left
plot(1e3*t_a,x(:,1),'linewidth',1.5)
ylabel('v_C (V)')
yyaxis right
plot(1e3*t_a,x(:,2),'linewidth',1.5)
xlabel('time (ms)')
ylabel('i_L (A)')
Warning: Imaginary parts of complex X and/or Y arguments ignored.
Warning: Imaginary parts of complex X and/or Y arguments ignored.

There’s something rotten in the state of Denmark. First of all, there are imaginary parts being ignored in the plots (see the warnings above), which is a bad sign (we could check them out to see how significant they are relative to the real parts). Even worse, as we can see in , a negative current \(i_L = i_D\) is predicted at certain times. This is literally impossible with a proper solution of this model because our diode model will not allow for negative current.

The most likely culprit here is the numerical solution. Our ODE is actually “stiff,” and would be better solved with a stiff solver like ode23s. Let’s try that.

[~,x] = ode23s( ...
    @(t,x) kafka_ode(t,x,u(t)), ... % ODE ... 
    t_a, ... % time array or span
    x0 ... % initial state
);

We now have our initial solution array x, which represents our state solution \(x(t)\). Let’s plot the solution.

figure
yyaxis left
plot(1e3*t_a,x(:,1),'linewidth',1.5)
ylabel('v_C (V)')
yyaxis right
plot(1e3*t_a,x(:,2),'linewidth',1.5)
xlabel('time (ms)')
ylabel('i_L (A)')

The result, shown in , had no “imaginary part” warnings and satisfies our knowledge that \(i_L(t) > 0\). And the response makes sense: the capacitor charges, the current stays non-negative, oscillations appear similar to the input frequency, and even three “cutoffs” of current when the system is tending toward negative current.

Part b: linear solution

Instead of using the usual linear simulation methods, for parallelism with the nonlinear solution, we begin by defining the linear state equation in the following function file (could be a local function).

type kafka_ode_lin.m
function dxdt = kafka_ode_lin(t,x,u)
    global R L C VTH IS
    A = [ ... 
      0, 1/C; ...
      -1/L, -(VTH/(2*IS)+R)/L ...
    ];
    B = [ ...
      0; 1/L ...
    ];
    dxdt = A*x + B*u; % linear state eq!
end

Out of curiousity, what are the eigenvalues?

A = [ ... 
      0, 1/C; ...
      -1/L, -(VTH/(2*IS)+R)/L ...
    ];
evals = eig(A);
disp(evals(1))
disp(evals(2))
     0

  -1.2500e+13

Interesting: a zero-eigenvalue and a large, negative, real one. The zero-eigenvalue corresponds to a capacitor voltage that grows without bound (in the model).

Now let’s use ode23 for the simulation.

options = odeset( ... % solver options
  'AbsTol',1e-13, ...
  'RelTol',1e-13 ...
);
[~,x_lin] = ode23s( ...
    @(t,x) kafka_ode_lin(t,x,u(t)), ... % ODE
    t_a, ... % time array
    x0, ... % initial state
    options ...
);

Now let’s remember that this is the solution for \(\bm{x}^*(t) = \bm{x}(t) - \bm{x}_o(t)\). Let’s “correct” for the operating point in the solution, even though it’s quite small.

xo = [0;IS]; % operating point
x_lin = x_lin + xo.';

We could plot the nonlinear and linear solutions together, but they’re so different we choose to plot the separately.

figure
yyaxis left
hold on
plot(1e3*t_a,x_lin(:,1),'linewidth',1.5)
xlabel('time (ms)')
ylabel('v_C (V)')
yyaxis right
plot(1e3*t_a,1e3*x_lin(:,2),'linewidth',1.5)
hold off
ylabel('i_L (mA)')
legend( ...
    'linear v_C', ...
    'linear i_L', ...
    'location','southeast' ...
)

part c: comparison

The plots of and show that the linear solution is absolutely terrible. This is to be expected with such strong nonlinearities and a state that varies significantly from the “operating point” assumed by the linearized model.

Problem 14.5 (HOOTENANNY)

Design a home rainwater catchment system and sprinkler distribution system. Most places, a surprising amount of water falls on a house’s roof throughout a year. Capturing it for irrigation can save water costs and reduce the environmental impact of watering lawns, plants, and gardens.

Design a home rainwater catchment and irrigation system. The design constraints are as follows.

  1. It should be designed for Olympia, Washington rainfall, as described in table 14.1.

  2. For a house, large tanks are unsightly. Instead, use a series of connected barrels.

After discussions with the customer, the following design requirements for the system are identified.

  1. It should be capable of distributing one inch of water per unit area June through September, even during drought conditions, during which there is half the average rainfall in the months March through September (see table 14.1.

  2. The roof area for collection is \(400\) square ft.

  3. The lawn area for distribution is \(600\) square ft.

  4. It should be low-maintenance.

  5. The distribution system should be capable of being “blown out” during winter months or it must be designed to handle sudden dips from \(33\) down to \(22\) deg F for up to two days.1

  6. When tanks are full, it should be able to gracefully dump excess water. If possible, designing it to refresh itself by dumping old water for new water is desired.

  7. It should be able to handle a heavy rain of \(1\) inch per hour via an overflow mechanism, but be able to handle a moderate rain of \(0.2\) inches per hour without requiring overflow (unless the tanks are full).

  8. It should be designed to be fed from typical house rain gutter downspouts.

  9. Distribution should be automated.

  10. Energy efficiency is desired. If possible, using tanks’ potential energy for distribution is desired. In this case, unconventional distribution networks are allowable (e.g. “drip” systems without conventional sprinkler heads that require high pressure). However, the distribution hardware should not be custom-designed.

  11. Commercially available parts are desired. Minimize the number of custom parts (zero is best).

The focus of the design problem is the sizing of the pipes, barrels, and mechanisms based on a dynamic system analysis.[^2]

It is highly recommended that you use the following Fourier Series fit to the Olympia drought rainfall data, presented as trigonometric series coefficient vectors a and b for easy definition in Matlab.2

w = 0.5236; % fundamental frequency
a0 = 3.579; % dc offset
a(1) = 4.144;
b(1) = 0.6244;
a(2) = 1.332;
b(2) = 0.07578;
a(3) = -0.07667;
b(3) = 0.03167;
a(4) = -0.2469;
b(4) = 0.0004836;
a(5) = -0.09448;
b(5) = 0.01735;
a(6) = 0.07417;
b(6) = -2.131e-06;
a(7) = -0.06748;
b(7) = -0.0124;
a(8) = -0.1235;
b(8) = -0.0002381;

A system model response to this input can be used to determine the system parameters, such as the number of barrels required. Do not forget to include the effect of distribution, which can be modeled as a negative source. Although we have the tools to perform the analysis analytically, it is highly recommended that a Matlab simulation is developed using ss to define the system and lsim to simulate the response. A frequency response analysis using bode may also prove useful. It may be possible to simply iteratively tweak design parameters until the simulation meets the requirements.

A thorough report is required. It is highly recommended that LaTeX is used. Thorough analysis, results, and design is required. All sizing and specific parts are required. Either an analytic or a numerical (simulation) demonstration of the design’s fulfillment of the requirements is required.

A good start on solving this is given in .

  1. A potential way to mitigate freezing is keeping the water in motion. Care must be taken not to create inadvertent ice skating rinks.↩︎

  2. The fit is an \(8\)-term Fourier series fit performed via Matlab’s fit function.↩︎

Online Resources for Section 14.6

No online resources.