... ... @@ -101,7 +101,7 @@ They describe (quasi-real) time series for wind power generation $$W(t)$$, solar \tilde{X}(\omega) = \int_0^T X(t) e^{\i \omega t} \,\ud t \, . \end{equation*} For all three regions, plot the energy spectrum $\left| \tilde{\Delta}(\omega) \right|^2$ as a function of $\left| \tilde{X}(\omega) \right|^2$ as a function of $\omega$. Discuss the relationship of these results with the findings obtained in (b)-(e). \item Normalize the time series to one, so that $$\expect{W} = \expect{S} = \expect{L} = 1$$. ... ...
 %% Cell type:markdown id: tags: # Energy System Modelling - Solutions to Tutorial II.1 # Energy System Modelling - Solutions to Tutorial I.2 %% Cell type:markdown id: tags: We use approximations to seasonal variations of wind and solar power generation $W(t)$ and $S(t)$ and load $L(t)$: $$W(t) = 1 + A_W \cos \omega t$$ $$S(t) = 1 - A_S \cos \omega t$$ $$L(t) = 1 + A_L \cos \omega t$$ The time series are normalized to $\langle{W}\rangle = \langle{S}\rangle = \langle{L}\rangle := \frac{1}{T} \int_0^T L(t) d t = 1$, and the constants have the values $$\omega = \frac{2\pi}{T}$$ $$T = 1 \text{ year}$$ $$A_W = 0.4$$ $$A_S = 0.75$$ $$A_L = 0.1$$ %% Cell type:markdown id: tags: ## Imports %% Cell type:code id: tags:  python import numpy as np import matplotlib.pyplot as plt import scipy.integrate as integrate from scipy.optimize import minimize %matplotlib inline  %% Cell type:markdown id: tags: ## Parameters %% Cell type:code id: tags:  python A_w = 0.4 A_s = 0.75 A_l = 0.1 T = 1 phi = 0  %% Cell type:markdown id: tags: ## Functions %% Cell type:code id: tags:  python def w(t, phi=0): return 1 + A_w * np.cos(2*t*np.pi/T-phi)  %% Cell type:code id: tags:  python def s(t): return 1 - A_s * np.cos(2*t*np.pi/T)  %% Cell type:code id: tags:  python def l(t): return 1 + A_l * np.cos(2*t*np.pi/T)  %% Cell type:code id: tags:  python def c(gamma): return 1 - gamma  %% Cell type:code id: tags:  python def mismatch_a(alpha, t, phi=phi): return alpha*w(t, phi) + (1-alpha)*s(t)-l(t)  %% Cell type:code id: tags:  python def mismatch_c(alpha, gamma, t): return gamma*(alpha*w(t) + (1-alpha)*s(t))+c(t)-l(t)  %% Cell type:code id: tags:  python def objective_a(alpha): return 1/T * integrate.quad(lambda t: mismatch_a(alpha,t,phi=phi)**2, 0, T)  %% Cell type:code id: tags:  python def objective_c(alpha): return 1/T * integrate.quad(lambda t: mismatch_c(alpha,gamma,t)**2, 0, T)  %% Cell type:markdown id: tags: ## Plots %% Cell type:code id: tags:  python x = np.arange(0,1,0.01)  %% Cell type:markdown id: tags: Availability time series %% Cell type:code id: tags:  python wind = plt.plot(x, w(x), label='wind') solar = plt.plot(x, s(x), label='solar') load = plt.plot(x, l(x), label='load') plt.legend() plt.show()  %% Output %% Cell type:markdown id: tags: Mismatch %% Cell type:code id: tags:  python alphas = [0, 0.25, 0.5, 0.75, 1] for a in alphas: plt.plot(x, mismatch_a(a,x),label='wind share = '+str(a)) plt.legend() plt.show()  %% Output %% Cell type:markdown id: tags: ## Solution %% Cell type:markdown id: tags: *** **(a) What is the seasonal optimal mix $\alpha$, which minimizes** $$\langle\left[ \alpha W(\cdot) + (1-\alpha) S(\cdot) - L(\cdot) \right]^2 \rangle = \frac1T \int_0^T \left[ \alpha W(t) + (1-\alpha) S(t) - L(t) \right]^2 \,\mathrm d t$$ > **Hint:** You can use the function [scipy.optimize.minimize](https://docs.scipy.org/doc/scipy-1.0.0/reference/generated/scipy.optimize.minimize.html) and use e.g. [method='nelder-mead'](https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method). %% Cell type:code id: tags:  python alpha = minimize(objective_a, x0=0, method='nelder-mead', options={'xtol': 1e-8, 'disp': True}).x alpha  %% Output Optimization terminated successfully. Current function value: 0.000000 Iterations: 37 Function evaluations: 74 0.73913043212890694 %% Cell type:markdown id: tags: *** **(b) How does the optimal mix change if we replace $A_L \to -A_L$?** %% Cell type:code id: tags:  python A_l *= (-1)  %% Cell type:code id: tags:  python alpha = minimize(objective_a, x0=0, method='nelder-mead', options={'xtol': 1e-8, 'disp': True}).x alpha  %% Output Optimization terminated successfully. Current function value: 0.000000 Iterations: 36 Function evaluations: 72 0.56521739196777387 %% Cell type:code id: tags:  python A_l *= (-1)  %% Cell type:markdown id: tags: *** **(c) Now assume that there is a seasonal shift in the wind signal $$W(t) = 1 + A_W \cos \left( \omega t - \phi \right).$$ Express the optimal mix $\alpha$ as a function of $\phi$.** > **Remark:** Note, that $\alpha\in [0,1]$ and you need to add this as bounds. > **Hint:** If you encounter problems (why?), try another optimisation algorithm [method='TNC'](https://en.wikipedia.org/wiki/Truncated_Newton_method) %% Cell type:code id: tags:  python phis = np.arange(0,2*np.pi,np.pi/100) alphas = [] for p in phis: phi = p bnds = [(0, 1)] alpha = minimize(objective_a, x0=0, method='TNC', bounds=bnds, options={'xtol': 1e-8}).x alphas.append(alpha) plt.plot(phis, alphas)  %% Output [] %% Cell type:markdown id: tags: *** **(d) A constant conventional power source $C(t) = 1 - \gamma$ is now introduced. The mismatch then becomes $$\Delta(t) = \gamma \left[ \alpha W(t) + (1-\alpha) S(t) \right] + C(t) - L(t)$$ Analogously to (a), find the optimal mix $\alpha$ as a function of $0 \leq \gamma \leq 1$, which minimizes $\langle{\Delta^2}\rangle$.** %% Cell type:code id: tags:  python gammas = np.arange(0,1,0.01) alphas = [] for g in gammas: gamma = g bnds = [(0, 1)] alpha = minimize(objective_c, x0=0, method='TNC', bounds=bnds, options={'xtol': 1e-8}).x alphas.append(alpha) plt.plot(gammas, alphas)  %% Output [] ... ...
 %% Cell type:markdown id: tags: # Energy System Modelling - Tutorial I SS 2018, Karlsruhe Institute of Technology, Institute for Automation and Applied Informatics *** %% Cell type:markdown id: tags: # Imports %% Cell type:code id: tags:  python import numpy as np import pandas as pd import matplotlib.pyplot as plt %matplotlib inline  %% Cell type:markdown id: tags: *** # Introductory Comments %% Cell type:markdown id: tags: ## Getting Help Executing cells with Shift-Enter and with h there is help. Help is available with . or load.sort_values() <- cursor between brackets, Shift- %% Cell type:markdown id: tags: ## Using one-dimensional arrays (Numpy and Pandas) %% Cell type:markdown id: tags: **Numpy** %% Cell type:code id: tags:  python a = np.arange(10) a  %% Output array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) %% Cell type:code id: tags:  python a[1:3]  %% Output array([1, 2]) %% Cell type:markdown id: tags: **Pandas** %% Cell type:code id: tags:  python s = pd.Series(np.random.random(3), index=['foo', 'bar', 'baz']) s  %% Output foo 0.748951 bar 0.286110 baz 0.274185 dtype: float64 %% Cell type:code id: tags:  python s["foo":"bar"]  %% Output foo 0.748951 bar 0.286110 dtype: float64 %% Cell type:markdown id: tags: ## Using two-dimensional arrays (Numpy and Pandas) %% Cell type:markdown id: tags: **Numpy** %% Cell type:code id: tags:  python np.random.random((3,5))  %% Output array([[0.28941884, 0.4004386 , 0.67028241, 0.40211146, 0.87918444], [0.42221308, 0.13808648, 0.08780219, 0.96793469, 0.62714239], [0.97390844, 0.83849364, 0.65386425, 0.15093854, 0.29176614]]) %% Cell type:markdown id: tags: **Pandas** %% Cell type:code id: tags:  python s = pd.DataFrame(np.random.random((3,5)), index=['foo', 'bar', 'baz']) s  %% Output 0 1 2 3 4 foo 0.390126 0.502056 0.958872 0.186866 0.860330 bar 0.042661 0.636768 0.233535 0.577435 0.909700 baz 0.957205 0.313162 0.348274 0.542685 0.698854 %% Cell type:code id: tags:  python s.mean()  %% Output 0 0.463331 1 0.483995 2 0.513560 3 0.435662 4 0.822961 dtype: float64 %% Cell type:markdown id: tags: *** # Problem I.1 The following data are made available to you on the __[coures homepage](https://nworbmot.org/courses/complex_renewable_energy_networks/)__: de_data.csv, gb_data.csv, eu_data.csv and alternatively wind.csv, solar.csv, load.csv They describe (quasi-real) time series for wind power generation $W(t)$, solar power generation $S(t)$ and load $L(t)$ in Great Britain (GB), Germany (DE) and Europe (EU). The time step is 1 h and the time series are several years long. > Remark: In this example notebook, we only look at Germany and the EU, Great Britain works in exactly the same way. %% Cell type:markdown id: tags: *** **Read Data** %% Cell type:code id: tags:  python de = pd.read_csv('tutorial_data/de_data.csv', parse_dates=True, index_col=0) eu = pd.read_csv('tutorial_data/eu_data.csv', parse_dates=True, index_col=0) gb = pd.read_csv('tutorial_data/gb_data.csv', parse_dates=True, index_col=0)  %% Cell type:code id: tags:  python wind = pd.read_csv('tutorial_data/wind.csv', parse_dates=True, index_col=0) solar = pd.read_csv('tutorial_data/solar.csv', parse_dates=True, index_col=0) load = pd.read_csv('tutorial_data/load.csv', parse_dates=True, index_col=0)  %% Cell type:markdown id: tags: Extra: Show the first 5 lines (header) of the German data: %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: Extra: Check that wind, solar and load files are just differently organized datasets and it's the same data: %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: *** **(a) Check that the wind and solar time series are normalized to ’per-unit of installed capacity’, and that the load time series is normalized to MW.** %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: *** **(b) Calculate the maximum, mean, and variance of the time series. ** %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: *** ** (c) For all three regions, plot the time series $W (t)$, $S(t)$, $L(t)$ for a winter month (January) and a summer month (July). ** %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: Extra: Also compare the wind between the different regions %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: *** **(d) For all three regions, plot the duration curve for $W(t)$, $S(t)$, $L(t)$.** > **Hint:** You might want to make use of the functions [.sort_values](https://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.sort_values.html) and [.reset_index](https://pandas.pydata.org/pandas-docs/version/0.23/generated/pandas.DataFrame.reset_index.html) > **Tip:** Go through the line de['wind'].sort_values(ascending=False).reset_index(drop=True).plot() dot by dot and note what happens to the output. %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: *** **(e) For all three regions, plot the probability density function of $W(t)$, $S(t)$, $L(t)$.** %% Cell type:markdown id: tags: There are two different methods: 1. [Histograms](https://en.wikipedia.org/wiki/Histogram) and 2. [Kernel density estimation (KDE)](https://en.wikipedia.org/wiki/Kernel_density_estimation). This [image](https://en.wikipedia.org/wiki/Kernel_density_estimation#/media/File:Comparison_of_1D_histogram_and_KDE.png) on the KDE page provides a good summary of the differences. You can do both with Panda! %% Cell type:markdown id: tags: First, let's look at the wind data: %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: Now, let's look at the solar data: %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: The solar data might be hard to see. Look at this in detail by limiting the y-axis to (0,2): %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: Finally, let's look at the load profile: %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: *** **(f) Apply a [(Fast) Fourier Transform](https://en.wikipedia.org/wiki/Fast_Fourier_transform) to the the three time series $X \in W(t), S(t), L(t)$:** $$\tilde{X}(\omega) = \int_0^T X(t) \;e^{i\omega t} \;\mathrm{d}t.$$ **For all three regions, plot the energy spectrum $\|\tilde{\Delta}(\omega)\|^2$ as a function of $\omega$. Discuss the relationship of these results with the findings obtained in (b)-(f).** %% Cell type:code id: tags:  python  **For all three regions, plot the energy spectrum $\|\tilde{X}(\omega)\|^2$ as a function of $\omega$. Discuss the relationship of these results with the findings obtained in (b)-(f).** %% Cell type:markdown id: tags: > **Hint:** To determine the frequencies rfffreq can be used, the argument d indicates the distance between two data points, 1h hour, which we specify as $\frac{1}{8760} a$, so that the frequencies come out in the unit $\frac{1}{a}$. > **Remark:** Use the function [numpy.fft.rfft](https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.rfft.html) and make sure you subtract the mean. > **Remark:** To determine the frequencies [numpy.fft.rfffreq](https://docs.scipy.org/doc/numpy-1.12.0/reference/generated/numpy.fft.rfftfreq.html) is used, the argument d indicates the distance between two data points, 1h hour, which we specify as $\frac{1}{8760} a$, so that the frequencies come out in the unit $\frac{1}{a}$. %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: *** **(g) Normalize the time series to one, so that $\langle{W}\rangle = \langle{S}\rangle = \langle{L}\rangle = 1$.** **Now, for all three regions, plot the mismatch time series** $$\Delta(t) = \gamma \alpha W(t) + \gamma (1 - \alpha) S(t) - L(t)$$ **for the same winter and summer months as in (c). Choose** 1. $\alpha \in \{0.0, 0.5, 0.75, 1.0\}$ with $\gamma = 1$, and 2. $\gamma \in \{0.5, 0.75, 1.0, 1.25, 1.5\}$ with $\alpha = 0.75$. **What is the interpretation of $\gamma$ and $\alpha$?** %% Cell type:markdown id: tags: Choose the country and alpha, gamma values and re-run: %% Cell type:code id: tags:  python d = de gamma = 1.0 alpha = 0  %% Cell type:markdown id: tags: Normalize the time series and calculate mismatch time series: %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: Plot the mismatch time series for the winter and summer months: %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: *** **(h) For all three regions, repeat (b)-(f) for the mismatch time series. What changed?** %% Cell type:markdown id: tags: **Data check** %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: **Time series plot** %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: **Duration curve** %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: **Probability density function** %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: **Fast Fourier Transform** %% Cell type:code id: tags:  python  ... ...
 %% Cell type:markdown id: tags: # Energy System Modelling - Tutorial II.1 # Energy System Modelling - Tutorial I.2 %% Cell type:markdown id: tags: We use approximations to seasonal variations of wind and solar power generation $W(t)$ and $S(t)$ and load $L(t)$: $$W(t) = 1 + A_W \cos \omega t$$ $$S(t) = 1 - A_S \cos \omega t$$ $$L(t) = 1 + A_L \cos \omega t$$ The time series are normalized to $\langle{W}\rangle = \langle{S}\rangle = \langle{L}\rangle := \frac{1}{T} \int_0^T L(t) d t = 1$, and the constants have the values $$\omega = \frac{2\pi}{T}$$ $$T = 1 \text{ year}$$ $$A_W = 0.4$$ $$A_S = 0.75$$ $$A_L = 0.1$$ %% Cell type:markdown id: tags: ## Parameters ## Imports %% Cell type:code id: tags:  python A_w = 0.4 A_s = 0.75 A_l = 0.1 T = 1 phi = 0 import numpy as np import matplotlib.pyplot as plt import scipy.integrate as integrate from scipy.optimize import minimize %matplotlib inline  %% Cell type:markdown id: tags: ## Imports ## Parameters %% Cell type:code id: tags:  python import numpy as np import matplotlib.pyplot as plt import scipy.integrate as integrate from scipy.optimize import minimize %matplotlib inline A_w = 0.4 A_s = 0.75 A_l = 0.1 T = 1 phi = 0  %% Cell type:markdown id: tags: ## Functions %% Cell type:code id: tags:  python def w(t, phi=0): return 1 + A_w * np.cos(2*t*np.pi/T-phi)  %% Cell type:code id: tags:  python def s(t): return 1 - A_s * np.cos(2*t*np.pi/T)  %% Cell type:code id: tags:  python def l(t): return 1 + A_l * np.cos(2*t*np.pi/T)  %% Cell type:code id: tags:  python def c(gamma): return 1 - gamma  %% Cell type:code id: tags:  python def mismatch(???): return ???  %% Cell type:code id: tags:  python def objective(alpha): return 1/T * integrate.quad(lambda t: mismatch(???)**2, 0, T)  %% Cell type:markdown id: tags: ## Plots %% Cell type:code id: tags:  python x = np.arange(0,1,0.01)  %% Cell type:markdown id: tags: Availability time series as shown on the worksheet. %% Cell type:code id: tags:  python wind = plt.plot(x, w(x), label='wind') solar = plt.plot(x, s(x), label='solar') load = plt.plot(x, l(x), label='load') plt.legend() plt.show()  %% Output %% Cell type:markdown id: tags: ## Problem I.2 %% Cell type:markdown id: tags: *** **(a) What is the seasonal optimal mix $\alpha$, which minimizes** $$\langle\left[ \alpha W(\cdot) + (1-\alpha) S(\cdot) - L(\cdot) \right]^2 \rangle = \frac1T \int_0^T \left[ \alpha W(t) + (1-\alpha) S(t) - L(t) \right]^2 \,\mathrm d t$$ > **Hint:** You can use the function [scipy.optimize.minimize](https://docs.scipy.org/doc/scipy-1.0.0/reference/generated/scipy.optimize.minimize.html) and use e.g. [method='nelder-mead'](https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method). %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: *** **(b) How does the optimal mix change if we replace $A_L \to -A_L$?** %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: *** **(c) Now assume that there is a seasonal shift in the wind signal $$W(t) = 1 + A_W \cos \left( \omega t - \phi \right).$$ Express the optimal mix $\alpha$ as a function of $\phi$.** > **Remark:** Note, that $\alpha\in [0,1]$ and you need to add this as bounds. > **Hint:** If you encounter problems (why?), try another optimisation algorithm [method='TNC'](https://en.wikipedia.org/wiki/Truncated_Newton_method) %% Cell type:code id: tags:  python  %% Cell type:markdown id: tags: *** **(d) A constant conventional power source $C(t) = 1 - \gamma$ is now introduced. The mismatch then becomes $$\Delta(t) = \gamma \left[ \alpha W(t) + (1-\alpha) S(t) \right] + C(t) - L(t)$$ Analogously to (a), find the optimal mix $\alpha$ as a function of $0 \leq \gamma \leq 1$, which minimizes $\langle{\Delta^2}\rangle$.** %% Cell type:code id: tags:  python  ... ...
 ... ... @@ -128,7 +128,7 @@ If you map the nodes to countries like \texttt{0=DK, 1=DE, 2=CH, 3=IT, 4=AT,5=CZ The linear power flow is given by \begin{equation} p_i = \sum_j \tilde{L}_{i,j}\theta_j \qquad \text{and} \qquad f_l = \frac{1}{x_l} \sum_i K_{i,l}\theta_i, \qquad \text{where} \qquad \tilde{L}_{i,j}= \sum_l = K_{i,l}\frac{1}{x_l} K_{j,l} p_i = \sum_j \tilde{L}_{i,j}\theta_j \qquad \text{and} \qquad f_l = \frac{1}{x_l} \sum_i K_{i,l}\theta_i, \qquad \text{where} \qquad \tilde{L}_{i,j}= \sum_l K_{i,l}\frac{1}{x_l} K_{j,l} \end{equation} is the weighted Laplacian. For simplicity, we assume identity reactance on all links $x_l = 1$. ... ...