The purpose of this exercise is to
- Derive a network residual model and implement it in MATLAB
- Solve network algebraic equations using MATLAB function fsolve
- Construct an DAE model of a network load in MATLAB
- Compute the time response using the implicit solver ode15i
STEP 1a
Download loadFlowOne.m, the residual function constructed in lecture and extend it to the full network. The first line of your function should be: [r,J] = dNetwork(x,s,z,p)
The inputs are defined as follows: %x: network voltages; x = [V2R V2I V3R V3I]; %s: network injections; s = [4 -3]; %z: dynamic parameters; z = [1 1]; %p: static parameters; p = [X1 X2 Vm Vp nsT Qm Qp nsC];
The static parameters have the following values: X1 = 0.03; X2 = .05; Vm = .98; Vp = 1.02; Qm = -.25; Qp = .25; nsT = 0.025; nsC = .2;
The outputs are defined as follows: %r: residual; %J: Jacobian of the residual; J = dr/dx
Next, derive the network residual model to compute the voltages at nodes 2 and 3. Note, there are 4 unknowns so you will need 4 independent residual equations. Follow the steps below.
- Apply the Kirchhoff current law to compute the power flowing out of node 2 through the two connected lines. Hint, you can eliminate V1 from the equation using the relation for a perfect transformer.
- The residual equations for node 2 are the real and imaginary parts of the power balance (sum of power flowing in and power flowing out). The power injection at this node coming from the turbine is P + iQ = s(1).
- Apply the Kirchhoff current law to compute the power flowing out of node 3 through the one connected line.
- The residual equations for node 3 are the real and imaginary parts of the power balance (sum of power flowing in and power flowing out). The power injection at this node coming from the load is P + iQ = s(2).
- Compute the Jacobians of the two residual equations with respect to the unknown variables.
- When you're finished your file should include a vector equation of the form r = f(x), where r has 4 entries. Then you should have a matrix equation of the form J = dr/dx(x), where J has dimension 4x4.
STEP 1b
Compute the voltages V2 and V3 using the following command: Vn = fsolve(@(x)dNetwork(x,s,z,p),[V2R0 V2I0 V3R0 V3I0],optimoptions('fsolve','Jacobian','on','tolx',1e-12,'tolfun',1e-6,'display','off'));
Use the flat start for the initial voltage values, i.e., V2R0 = 1; V3R0 = 1; V2I0 = 0; V3I0 = 0;
The result should be: Vn = [0.9846 0.0300 0.9651 -0.1229];
STEP 2a
Download odefunBall.m, the residual function defining the DAE system for a pendulum. Use this as a template to construct a function for the network load. Notice, the model of the load your received in lecture included only 3 residues but the state below has five components. Hence you will need to add two fictitious residues for states 5 and 6. A common assumption is that the states not described in the module remain constant, in other words,
r(4) = d(4)
r(5) = d(5);
The first line of you function should be:
r = dLoad(t,x,d,z,p)
where % t = time % x: the continuous state vector; x = [xd Pd Qd V3R V3I]; %z: discrete state vector; z = []; %p: parameter vector; p = [Tp Ps Pt]; where Tp = 5; Ps = 3; Pt = .4;
STEP 2b
Compute the time response of the load using the ode15i solver. First you will need to call the function decic to compute the initial state:
[x0,d0] = decic(@(t,x,d)dLoad(t,x,d,z,p),t0,[xd0 Pd0 Qd0 V3R0 V3I0],[0 0 0 1 1],[0 0 0 0 0],[0 0 0 0 0]); where t0 = 0; xd0 = 2.6214; Pd0 = 3; Qd0 = 0; and the real and imaginary parts of voltage at node 3 are taken from Part1 of this exercise.
Next, compute the the response as follows:
sol = ode15i(@(t,x,d)odefunBall(t,x,d,z,p),[0 15],x0p,d0); where d0 is taken from above and x0p is obtained from x0 as follows;
x0p = x0; x0p(4:5) = x0p(4:5)*.9;
This represents a 10% change in the voltage at node 3.
The time response should match the figure below:
