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

add simulation example for lecture on uncertain system analysis

parent 7fb72e44
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Index.ipynb
+ 8
− 0
View file @ f21a0c44
%% Cell type:markdown id: tags:
# Modeling and Simulation of Complex Power Systems
## Simulation Examples
On this platform you find the simulation examples based on Jupyter Notebooks. They are presented in the course **Modeling and Simulation of Complex Power Systems** by **Prof. Monti** at the **RWTH Aachen University**. Reviewing and modifying the simulation examples shall support the students in understanding the lecture content.
The following notebooks are currently available:
- [Lecture 2 - Modified Nodal Analysis](./lectures/02_NA_MNA/VS_CS_R4.ipynb)
- [Lecture 3 - Resistive Companion](./lectures/03_ResistiveCompanion/VS_R2L3.ipynb)
- [Lecture 4 - Nonlinear Resistive Companion](./lectures/04_NLResistiveCompanion/NL_RC.ipynb)
- [Lecture 5 - State Space Equations](./lectures/05_StateSpace/StateEq_ASMG.ipynb)
- [Lecture 8 - Dynamic Phasors](./lectures/08_DecoupledELMESim/VS_RL1.ipynb)
- [Lecture 9 - Diakoptics](./lectures/09_Diakoptics/Diakoptics.ipynb)
- [Lecture 9 - LIM](./lectures/09_LIM/LIM.ipynb)
- [Lecture 11 - Uncertain System Analysis](./lectures/11_UncertainSystemAnalysis/UncertainSystemAnalysis.ipynb)
## Helping Material
### Notebooks
Learn how to work with the notebooks by accessing the extensive information under **Help** in the navigation bar.
As a first starting point you might refer to **Help>Notebook Reference**, where you can find **Notebook Examples** in the user documentation. Important basics are covered in the following sections:
- Notebook Basics
- Running Code
- Markdown Cells
### Python
If you need an introduction to the Python language, the **Beginners Guide** under **Help>Python Reference** links to lots of tutorials depending on your previous knowledge.
### IPython
As a more advanced user, you can benefit from the fact that Jupyter Notebooks run an **IPython** kernel, which comes with additional features with respect to the standard **Python** kernel. To learn more about these features, you can follow under **Help>IPython Reference** the link **IPython documentation** and have a look into **Tutorial>Introducing IPython**.
### Contact
Any questions can be adressed to
[acs-teaching-mscps@eonerc.rwth-aachen.de](mailto:acs-teaching-mscps@eonerc.rwth-aachen.de)
%% Cell type:code id: tags:
``` python
```
......
lectures/11_UncertainSystemAnalysis/Circuit_RL.png 0 → 100644
+ 0
− 0
View file @ f21a0c44
lectures/11_UncertainSystemAnalysis/Circuit_RL.png

76.3 KiB

lectures/11_UncertainSystemAnalysis/UncertainSystemAnalysis.html 0 → 100644
+ 12196
− 0
View file @ f21a0c44
Source diff could not be displayed: it is too large. Options to address this: view the blob.
lectures/11_UncertainSystemAnalysis/UncertainSystemAnalysis.ipynb 0 → 100644
+ 393
− 0
View file @ f21a0c44
%% Cell type:markdown id: tags:
# MSP Simulation Example
# Uncertain System Analysis - Classical Polynomial Chaos
%% Cell type:markdown id: tags:
## Sample Circuit
%% Cell type:markdown id: tags:
<img src="Circuit_RL.png" width="300" align="left">
%% Cell type:markdown id: tags:
$L$: $5 mH$
$E$: $2 V$
Uncertain parameter in the system is the resistance with the uniform distribution:
$R$: $R_0 \pm \Delta R $
$R_0$: $0.4 \Omega$
$\Delta R$: $0.05 \Omega$
%% Cell type:markdown id: tags:
### Deterministic Solution of the Circuit
Euler-Forward integration method:
$I(k) = (1-T_s*R/L) * I(k-1) + Ts/L*E$
%% Cell type:markdown id: tags:
## Circuit and Simulation Setup
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(sign=' ')
# Circuit parameters:
E = 2.0
L = 5.0e-3
R0 = 0.4
R1 = 0.05
# Total simulation time
T_total = 0.08
# Simulaiton time step
Ts = 0.1e-3
# Number of simulation time steps
npoint = int(np.round(T_total/Ts))
print('Total simulation time: ' + str(T_total))
print('\nSimulation time step: ' + str(Ts))
```
%% Output
Total simulation time: 0.08
Simulation time step: 0.0001
%% Cell type:markdown id: tags:
## Monte Carlo Simulation Method
%% Cell type:code id: tags:
``` python
# Number of samples for Monte Carlo simultion runs
N = 40
# Generate random variable
rand_var = np.random.uniform(-1, 1, N)
print('Random variable: \n' + str(rand_var))
# R is a random parameter with the uniform distribution R= R0 + R1*rand_var
R = np.array([ R0 + R1*rand_var ])
print('\n Resistor R: \n' + str(R))
# Solution vector for current I
# Each row refers to a solution vector obtained in one Monte Carlo simulation run
I = np.zeros((N, npoint))
# Monte Carlo simulation loop (for each sample of random parameter R)
for j in np.arange(0,N-1):
# Value of random parameter R for this Monte Carlo simulation run
Rj = R[0, j]
# Time loop
# Initial condition
I[j,0]=0;
for i in np.arange(1,npoint):
# Euler Forward used for discretization
I[j,i] = (1-Ts*Rj/L) * I[j,i-1] + Ts/L*E
```
%% Output
Random variable:
[-8.59029337e-01 -4.60466494e-01 5.91887395e-02 -1.46515604e-01
-7.15203525e-02 -8.06258266e-01 -3.31665853e-01 6.65056465e-01
-9.24720924e-01 -9.64876536e-01 9.88020137e-01 1.59075805e-01
-4.12507969e-01 7.24957527e-01 -8.55870483e-01 8.05672962e-01
-2.77408910e-01 1.22389746e-01 9.36528164e-01 -2.44663481e-01
-4.30228245e-01 2.47867440e-01 -1.88870745e-04 6.75921027e-01
1.48082059e-01 -1.05219493e-01 -3.15744588e-01 -4.68967536e-01
9.69742135e-01 4.99619187e-01 -6.44014720e-01 1.72531036e-01
-5.20327536e-01 -7.23361683e-01 -3.75604850e-02 -6.59856990e-01
-7.60130822e-01 4.57536217e-01 4.29482901e-01 2.68400232e-01]
Resistor R:
[[ 0.35704853 0.37697668 0.40295944 0.39267422 0.39642398 0.35968709
0.38341671 0.43325282 0.35376395 0.35175617 0.44940101 0.40795379
0.3793746 0.43624788 0.35720648 0.44028365 0.38612955 0.40611949
0.44682641 0.38776683 0.37848859 0.41239337 0.39999056 0.43379605
0.4074041 0.39473903 0.38421277 0.37655162 0.44848711 0.42498096
0.36779926 0.40862655 0.37398362 0.36383192 0.39812198 0.36700715
0.36199346 0.42287681 0.42147415 0.41342001]]
%% Cell type:markdown id: tags:
### Plots
%% Cell type:code id: tags:
``` python
# Plots
# Time vector
t = np.arange(0, npoint)*Ts
plt.figure(figsize=(8,6))
for j in np.arange(0,N-1):
plt.plot(t,I[j,:], linewidth=1)
plt.xlabel('Time [s]')
plt.ylabel('Current [A]')
plt.title('Monte Carlo simulation results')
plt.show()
```
%% Output
%% Cell type:markdown id: tags:
## Classical Polynomial Chaos
%% Cell type:markdown id: tags:
### Polynomial Chaos (PC) expansion
PC expansion of $R$:
$R = \sum_{i=0}^{P}R_i\Phi_i(\xi) = R_0\Phi_0(\xi) + R_1\Phi_1(\xi)$
PC expansion of $i(t)$:
$I = \sum_{i=0}^{P}I_i(t)\Phi_i(\xi) = I_0(t)\Phi_0(\xi) + I_1(t)\Phi_1(\xi)$
Circuit solution:
$I(k) = (1-T_s*R/L) * I(k-1) + Ts/L*E = I(k-1) - T_s/L*R*I(k-1) + Ts/L*E$
Replacing $R$ and $I$ to the solution equation:
$\sum_{i=0}^{P}I_i(k)\Phi_i = \sum_{i=0}^{P}I_i(k-1)\Phi_i - T_s/L\sum_{i=0}^{P}\sum_{j=0}^{P}R_i*I_j(k-1)\Phi_i\Phi_j + T_s/L*E$
Applying Galerkin projection on the PC basis and by replcing integrals with inner products:
$\sum_{i=0}^{P}I_i(k)<\Phi_i\Phi_s> = \sum_{i=0}^{P}I_i(k-1)<\Phi_i\Phi_s> - T_s/L\sum_{i=0}^{P}\sum_{j=0}^{P}R_i*I_j(k-1)<\Phi_i
\Phi_j\Phi_s> + T_s/L*E$
Inner product of two orthogonal polynomials can be replaced by the following identity:
$<\Phi_i\Phi_s> = <\Phi_s^2>\delta_{is}$
where $\delta_{is}$ is Kronecker delta.
The following can be obtained:
$\sum_{i=0}^{P}I_i(k)<\Phi_s^2>\delta_{is} = \sum_{i=0}^{P}I_i(k-1)<\Phi_s^2>\delta_{is} - T_s/L\sum_{i=0}^{P}\sum_{j=0}^{P}R_i*I_j(k-1)<\Phi_i\Phi_j\Phi_s> + T_s/L*E$
Kronecker delta reduces elements in summations in the following:
$I_s(k)<\Phi_s^2> = I_s(k-1)<\Phi_s^2> - T_s/L\sum_{i=0}^{P}\sum_{j=0}^{P}R_i*I_j(k-1)<\Phi_i\Phi_j\Phi_s> + T_s/L*E$
The PC expansion of the circuit solution is:
$I_s(k) = I_s(k-1) - T_s/L\frac{1}{<\Phi_s^2>}\sum_{i=0}^{P}\sum_{j=0}^{P}R_i*I_j(k-1)<\Phi_i\Phi_j\Phi_s> + T_s/L*E$
%% Cell type:markdown id: tags:
#### Matching of PDF types and orthogonal polynomials
Uniform distribution of parameter -> Legendre polynomial
Order of PC expansion:
$P = 1$
Non-zero inner products for Legendre polynomials:
$<\Psi_0\Psi_0> = 1$
$<\Psi_1\Psi_1> = 1/3$
$<\Psi_0\Psi_0\Psi_0> = 1$
$<\Psi_0\Psi_1\Psi_1> = <\Psi_1\Psi_0\Psi_1> = <\Psi_1\Psi_1\Psi_0> = 1/3$
The first order PC expansion of the circuit solution is the following:
$s=0$
$I_0(k) = (1-T_s*R_0/L) * I_0(k-1) - T_s/L/3*R_1*I_1(k-1) + Ts/L*E$
$s=1$
$I_1(k) = -T_s/L*R_1*I_0(k-1) - (1-T_s*R_0/L)*I_1(k-1)$
%% Cell type:markdown id: tags:
#### Calculation of PC coefficients
%% Cell type:code id: tags:
``` python
# P-th order PC expansioin
P=1
# Number of coefficients in PC expansioin
M=P+1
# Solution vector for PC
# Each row refers to a solution vector obtained for a PC coefficient
I_pct = np.zeros((M, npoint))
# Matrix for solution equation
Gh = np.array([ [(1-Ts*R0/L), -Ts*R1/L/3],
[-Ts*R1/L, (1-Ts*R0/L)]])
# Simulation time loop
for k in np.arange(1,npoint):
# Euler Forward used for discretization
I_pct_k = np.reshape(np.matmul(Gh, I_pct[:,k-1]), (M,1)) + np.array([ [Ts/L*E], [0]])
I_pct[0,k] = I_pct_k[0,0]
I_pct[1,k] = I_pct_k[1,0]
```
%% Cell type:markdown id: tags:
### Reconstruction of uncertain variable based on PC coefficients
%% Cell type:markdown id: tags:
Legendre polynomial of order $1$:
$\Phi(\xi) = \Phi_0(\xi) + \Phi_1(\xi) = 1 + \xi$
$I = \sum_{i=0}^{P}I_i\Phi_i(\xi) = I_0\Phi_0(\xi) + I_1\Phi_1(\xi) = I_0 + I_1\xi$
%% Cell type:markdown id: tags:
### Plots
%% Cell type:code id: tags:
``` python
# Time vector
t = np.arange(0, npoint)*Ts
# Solution vector
i_pct = np.zeros((M, npoint))
# Number of reconstructions
K = 15
plt.figure(figsize=(8,6))
for j in np.arange(0,K-1):
# Reconstruction of random behaviour of I based on PC expansioin coefficients
i_pct = I_pct[0,:] + I_pct[1,:]*rand_var[j]
plt.plot(t,i_pct, linewidth=1)
plt.xlabel('Time [s]')
plt.ylabel('Current [A]')
plt.title('Polynomial Chaos simulation results')
plt.show()
```
%% Output
%% 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