Skip to content
Snippets Groups Projects
Commit 1b0d9ec3 authored by Jan Dinkelbach's avatar Jan Dinkelbach
Browse files

improve plotting in case 2

parent 240e776d
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
SLEW_Case2_Electromechanical_Phenomena.ipynb
+ 17
− 41
View file @ 1b0d9ec3
%% Cell type:markdown id:cf1f4a01-23c9-483d-9bf3-0a80ed440ae0 tags:
%% Cell type:markdown id: tags:
<div>
<img src="figures/slew_logo.png" style="float: right;height: 15em;">
</div>
# SLEW Case 2 - Electromechanical Dynamics: Small-signal Stability
## Theoretical background
### The Linearized Swing Equation
The swing equation is the fundamental equation that determines the rotor dynamics of a synchoronous machine.
$$M\frac{d^{2}\delta}{dt^{2}} = P_{m} - P_{e}(\delta) - D\frac{d\delta}{dt}$$
For small signal stability analysis, the swing equation can be linearized in the vicinity of the operating point $\delta = \delta_{0}$:
$$M\frac{d^{2}\Delta\delta}{dt^{2}} + D\frac{d\Delta\delta}{dt} + K_{e^{'}}\Delta\delta = 0 $$
where $K_{e^{'}}$ is the transient synchoronizing power coefficient at $\delta = \delta_{0}$:
$$K_{e^{'}} = \frac{\partial P_{e^{'}}}{\partial \delta^{'}}\mid_{\delta = \delta_{0}}$$
The characteristic equation and the roots of the linearized swing equation are as follows:
$$\lambda^{2} + \frac{D}{M}\lambda + \frac{K_{e^{'}}}{M} = 0 $$
$$\lambda_{1,2} = -\frac{D}{2M} \pm \sqrt{\left(\frac{D}{2M}\right)^{2} - \frac{K_{e^{'}}}{M}} $$
%% Cell type:markdown id:a620deac-bea3-4f34-841a-b9d7e41fe1ae tags:
%% Cell type:markdown id: tags:
***
## Case description
%% Cell type:markdown id:40672b8b-585f-4448-9dac-523d4f33287c tags:
%% Cell type:markdown id: tags:
A $200 MVA$ round-rotor generator is connected to an infinite bus system of nominal frequency of $50 Hz$ through a step-up transformer with reactance of $0.13 p.u.$ and a transmission line with reactance of $0.17 p.u.$ . The generator’s transient reactance is $x'_{d} = 0.23 p.u.$ and the inertia constant is $H = 4 s$. The synchoronous generator mechanical power is set to $0.6 p.u.$ and the steady-state emf $E = 1.1 p.u.$
%% Cell type:markdown id:ab2be9d0-04eb-47cd-988d-3af3cde391bd tags:
%% Cell type:markdown id: tags:
***
## Case tasks
%% Cell type:markdown id:4947469f-3514-4edc-a9d9-4d3157f36513 tags:
%% Cell type:markdown id: tags:
1. If at the given operating point the system eigen-values are as follows:
$$\lambda_{1,2} = -0.06254 \pm 8.8318j $$
a.) Determine whether the system response following a small disturbance is overdamped, critically damped or underdamped.<br>
b.) Verify you findings by calculating the damping coefficient of the synchronous generator and then applying the damping value in VILLASweb and generating the corresponding power and angle plots.<br>
<br>
2. If the system response following a small disturbance is critically damped. <br>
a.) Determine the synchoronous generator damping coefficient.<br>
b.) Apply the calculated damping value in VILLASweb and generate the corresponding power and angle plots.<br>
c.) Increase the damping coefficient value and observe the impact of such increase.<br>
%% Cell type:markdown id:68a26943-3791-4b08-a38e-a4a1454d24b7 tags:
%% Cell type:markdown id: tags:
**You can process the results from VILLASweb (as required in 1 and 2) using the prepared notebook cells below.**
%% Cell type:markdown id:8b454ae8-71a1-4017-9c5a-cf3efc494e40 tags:
%% Cell type:markdown id: tags:
***
## Case solutions
%% Cell type:markdown id:fa5f07a4-23a0-4846-8fc2-d9a3eeae3b8a tags:
%% Cell type:markdown id: tags:
#### First Setup your Python for post-processing
%% Cell type:markdown id:032049e7-1108-4fb1-9057-3b6c522a81bd tags:
%% Cell type:markdown id: tags:
Import relevant Python packages by executing the following cell:
%% Cell type:code id:fa7d426a-0908-454d-8b0d-0dc64eb785ca tags:
%% Cell type:code id: tags:
``` python
from IPython.display import display
from villas.web.result import Result
from pprint import pprint
import tempfile
import os
from villas.dataprocessing.readtools import *
from villas.dataprocessing.timeseries import *
import matplotlib.pyplot as plt
import scipy.io as sio
outputs = sio.loadmat('outputs/case2_output.mat', simplify_cells=True)
%matplotlib widget
```
%% Cell type:markdown id:9be353cf-624a-4e88-a427-f94091c651bb tags:
%% Cell type:markdown id: tags:
### Subtask 1.b
To check your results, enter your solution for $D$ rounded to 4 digits behind the comma:
%% Cell type:code id:cd8c51d0-ca49-426b-8b86-af774b6717b7 tags:
%% Cell type:code id: tags:
``` python
D=input('D in Subtask 1.b:')
print('Your result is', 'correct!' if round(float(D),4)==round(outputs['D_T1b'],4) else 'incorrect. You should double-check your calculations.')
```
%% Cell type:markdown id:ab551c32-ef6f-47d3-b839-57e2ce506863 tags:
%% Cell type:markdown id: tags:
#### Download simulation results from VILLASweb
%% Cell type:markdown id:e70b8ccc-8d97-4827-ba45-251be9c6ac82 tags:
%% Cell type:markdown id: tags:
Enter VILLASweb under https://slew.rwth-aachen.de/.
Adjust the parameter $D$ and run a corresponding simulation.
Then, enter your code snippet for result download in the following cell:
%% Cell type:code id:5c6f4d63-2c1c-49aa-acde-073f1cffd7d5 tags:
%% Cell type:code id: tags:
``` python
# Access result using token
r = Result(1, 'xyz', endpoint='https://slew.rwth-aachen.de')
```
%% Cell type:markdown id:85d62e17-9df8-421c-8b99-aac7be09bf09 tags:
%% Cell type:markdown id: tags:
#### Plot the results
%% Cell type:code id:ec357284-6407-41bb-9d54-3ca6c0532fa9 tags:
%% Cell type:code id: tags:
``` python
results_file_name='SP_SynGenTrStab_SMIB_Fault_SlewCase2_JsonSyngenParams_SP.csv'
# Get file by name
f = r.get_file_by_name(results_file_name)
# Create temp dir
tmpdir = tempfile.mkdtemp()
results_file_path = os.path.join(tmpdir,results_file_name)
# Load files as bytes
f = f.download(dest=results_file_path)
# Get time series object from csv
ts_res1 = read_timeseries_csv(results_file_path, print_status=True)
# Plot the currents
plt.figure(figsize=(8,12))
name_list = ['Rotor Angle','Output Power']
plt.subplot(2,1,1)
plt.ylabel("Angle (degree)")
plt.xlabel("time (s)")
plt.plot(ts_res1['delta_r_gen'].time, ts_res1['delta_r_gen'].values, label='Rotor Angle task 1.b', color='C0')
plt.xlim([1, 30])
plt.legend(loc='upper right')
plt.subplot(2,1,2)
plt.ylabel("Power (MW)")
plt.xlabel("time (s)")
plt.plot(ts_res1['P_elec'].time, ts_res1['P_elec'].values/1e6, label='Output Power task 1.b', color='C0')
plt.xlim([1, 30])
plt.legend(loc='upper right')
plt.show()
```
%% Cell type:markdown id:87835d8e-afcd-4a23-8b36-596eaaca7b92 tags:
%% Cell type:markdown id: tags:
### Subtask 2.a
To check your results, enter your solution for $D$ rounded to 4 digits behind the comma:
%% Cell type:code id:56d9f0db-b158-4fc3-a9b0-5a2d8931da9b tags:
%% Cell type:code id: tags:
``` python
D=input('D in Subtask 2.a:')
print('Your result is', 'correct!' if round(float(D),4)==round(outputs['D_T2a'],4) else 'incorrect. You should double-check your calculations.')
```
%% Cell type:markdown id:0848fc38-50ad-421a-8442-9f8b9c26a1e0 tags:
%% Cell type:markdown id: tags:
### Subtask 2.b
%% Cell type:markdown id:3c0fd278-b58a-4af7-a101-4a52e1a32c94 tags:
%% Cell type:markdown id: tags:
#### Download simulation results from VILLASweb
%% Cell type:markdown id:af0ec7ca-9cf7-4986-9ae9-ab762a792b5f tags:
%% Cell type:markdown id: tags:
Enter VILLASweb under https://slew.rwth-aachen.de/.
Adjust the parameter $D$ and run a corresponding simulation.
Then, enter your code snippet for result download in the following cell:
%% Cell type:code id:8331e105-8b66-4479-9caf-6ec038c1bce0 tags:
%% Cell type:code id: tags:
``` python
# Access result using token
r = Result(1, 'xyz', endpoint='https://slew.rwth-aachen.de')
```
%% Cell type:markdown id:f9879ce7-72ae-4306-9d65-8f8604484c43 tags:
%% Cell type:markdown id: tags:
#### Plot the results
%% Cell type:code id:ccad2758-d711-46c9-a54b-a788efb563eb tags:
%% Cell type:code id: tags:
``` python
results_file_name='SP_SynGenTrStab_SMIB_Fault_SlewCase2_JsonSyngenParams_SP.csv'
# Get file by name
f = r.get_file_by_name(results_file_name)
# Create temp dir
tmpdir = tempfile.mkdtemp()
results_file_path = os.path.join(tmpdir,results_file_name)
# Load files as bytes
f = f.download(dest=results_file_path)
# Get time series object from csv
ts_res2 = read_timeseries_csv(results_file_path, print_status=False)
# Plot the currents
plt.figure(figsize=(8,12))
name_list = ['Rotor Angle','Output Power']
plt.subplot(2,1,1)
plt.ylabel("Angle (degree)")
plt.xlabel("time (s)")
plt.plot(ts_res1['delta_r_gen'].time, ts_res1['delta_r_gen'].values, label='Rotor Angle task 1.b', color='C0')
plt.plot(ts_res2['delta_r_gen'].time, ts_res2['delta_r_gen'].values, label='Rotor Angle task 2.b', color='C1', linestyle=':')
plt.plot(ts_res2['delta_r_gen'].time, ts_res2['delta_r_gen'].values, label='Rotor Angle task 2.b', color='C1')
plt.plot(ts_res1['delta_r_gen'].time, ts_res1['delta_r_gen'].values, label='Rotor Angle task 1.b', color='C0', linestyle=':')
plt.xlim([1, 30])
plt.legend(loc='upper right')
plt.subplot(2,1,2)
plt.ylabel("Power (MW)")
plt.xlabel("time (s)")
plt.plot(ts_res1['P_elec'].time, ts_res1['P_elec'].values/1e6, label='Output Power task 1.b', color='C0')
plt.plot(ts_res2['P_elec'].time, ts_res2['P_elec'].values/1e6, label='Output Power task 2.b', color='C1', linestyle=':')
plt.plot(ts_res2['P_elec'].time, ts_res2['P_elec'].values/1e6, label='Output Power task 2.b', color='C1')
plt.plot(ts_res1['P_elec'].time, ts_res1['P_elec'].values/1e6, label='Output Power task 1.b', color='C0', linestyle=':')
plt.xlim([1, 30])
plt.legend(loc='upper right')
plt.show()
```
%% Cell type:markdown id:8e4f8736-d4a5-4d5b-8d56-fee4b7afb499 tags:
%% Cell type:markdown id: tags:
### Subtask 2.c
%% Cell type:markdown id:88f14679-151b-47d3-bf9e-0c912427e504 tags:
%% Cell type:markdown id: tags:
#### Download simulation results from VILLASweb
%% Cell type:markdown id:0740d062-5e8f-4eb9-88bc-a6eea0003391 tags:
%% Cell type:markdown id: tags:
Enter VILLASweb under https://slew.rwth-aachen.de/.
Adjust the parameter $D$ and run a corresponding simulation.
Then, enter your code snippet for result download in the following cell:
%% Cell type:code id:e490e568-1db0-4fb6-942a-30fafecee4bd tags:
%% Cell type:code id: tags:
``` python
# Access result using token
r = Result(1, 'xyz', endpoint='https://slew.rwth-aachen.de')
```
%% Cell type:markdown id:40ee9509-a952-42d8-a67d-4901ba2d95bd tags:
%% Cell type:markdown id: tags:
#### Plot the results
%% Cell type:code id:36505e96-7b81-41df-847c-a024269beb05 tags:
%% Cell type:code id: tags:
``` python
results_file_name='SP_SynGenTrStab_SMIB_Fault_SlewCase2_JsonSyngenParams_SP.csv'
# Get file by name
f = r.get_file_by_name(results_file_name)
# Create temp dir
tmpdir = tempfile.mkdtemp()
results_file_path = os.path.join(tmpdir,results_file_name)
# Load files as bytes
f = f.download(dest=results_file_path)
# Get time series object from csv
ts_res3 = read_timeseries_csv(results_file_path, print_status=False)
# Plot the currents
plt.figure(figsize=(8,12))
name_list = ['Rotor Angle','Output Power']
plt.subplot(2,1,1)
plt.ylabel("Angle (degree)")
plt.xlabel("time (s)")
plt.plot(ts_res2['delta_r_gen'].time, ts_res2['delta_r_gen'].values, label='Rotor Angle task 2.b', color='C1')
plt.plot(ts_res3['delta_r_gen'].time, ts_res3['delta_r_gen'].values, label='Rotor Angle task 2.c', color='C2', linestyle=':')
plt.plot(ts_res3['delta_r_gen'].time, ts_res3['delta_r_gen'].values, label='Rotor Angle task 2.c', color='C2')
plt.plot(ts_res2['delta_r_gen'].time, ts_res2['delta_r_gen'].values, label='Rotor Angle task 2.b', color='C1', linestyle=':')
plt.xlim([1, 30])
plt.legend(loc='upper right')
plt.subplot(2,1,2)
plt.ylabel("Power (MW)")
plt.xlabel("time (s)")
plt.plot(ts_res2['P_elec'].time, ts_res2['P_elec'].values/1e6, label='Output Power task 2.b', color='C1')
plt.plot(ts_res3['P_elec'].time, ts_res3['P_elec'].values/1e6, label='Output Power task 2.c', color='C2', linestyle=':')
plt.plot(ts_res3['P_elec'].time, ts_res3['P_elec'].values/1e6, label='Output Power task 2.c', color='C2')
plt.plot(ts_res2['P_elec'].time, ts_res2['P_elec'].values/1e6, label='Output Power task 2.b', color='C1', linestyle=':')
plt.xlim([1, 30])
plt.legend(loc='upper right')
plt.show()
```
%% Cell type:code id: tags:
``` python
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment