diff git a/.gitignore b/.gitignore
index ac3067842fab78d5e2ba485a3977f781f1bbf92f..0fa8a8b83f0d27695f4f0f49866e1585898a9eca 100644
 a/.gitignore
+++ b/.gitignore
@@ 1,8 +1,9 @@
# Ignore documentation
doc/Doxygen
doc/build
+src/
n# Ignore data and results
+# Ignore data and results
results/
# Ignore irrelevant personal files in
diff git a/.gitlabci.yml b/.gitlabci.yml
index 85a02f48596a70dfb2a8d6430a10d6155f807666..915a96fba0c7d412ee2fd523152efee495e30742 100644
 a/.gitlabci.yml
+++ b/.gitlabci.yml
@@ 55,7 +55,7 @@ pages:
script:
# Install texlive
 aptget y update
  aptget y install texlivefull
+ #  aptget y install texlivefull
# Install setuptools and sphinx
 pip install setuptools sphinx sphinx_rtd_theme
 pip install e git+https://github.com/bitprophet/releases/#egg=releases
diff git a/doc/source/user_guide/demand.rst b/doc/source/user_guide/demand.rst
index 8b9b78feb203a3ca8c9ad7b7a68faf8f8a1440f4..8ff57a7519a0696e89d246f533dfee73bdd1e2cf 100644
 a/doc/source/user_guide/demand.rst
+++ b/doc/source/user_guide/demand.rst
@@ 16,16 +16,6 @@ we adopt a Bayesian approach (O'Hagan, 1994) allowing one not only
to predict the mean of the demand conditioned on the temperature,
but also the conditional distribution around this mean.
Let the electricity demand :math:`D_i` for the region :math:`i`
be given by

.. math::
 D_i(T_i) = f_i(T_i) + \epsilon.
 :label: demandModel

where :math:`f_i, 1 \le i \le N` is some realvalued function
of the temperature :math:`T_i` in the region :math:`i`
and the residue :math:`\epsilon` is also to be modeled.
Model design
============
@@ 35,18 +25,7 @@ Electric heating is switched on only for lower temperatures,
while air conditioning is switched on only for higher temperatures.
Once an individual appliance is switched on, its electricity
consumption is to a first approximation linear in the temperature.
Assuming, that all consumers behave in the same way and that
a consumer switches the heater (air conditioning) for a
constant temperature threshold :math:`T_H` (:math:`T_C`),
we are led to define the functions :math:`f_i` as the
following piecewiselinear function
.. math::
 f_i(T_i) = w_H \Theta(T_H  T_i) (T_H  T_i)
 + w_C \Theta(T_H  T_i) (T_i  T_H) + w_0,
 :label: demandFunction

where :math:`\Theta` is the Heaviside distribution.
In fact, the relationship between the demand and the ambient temperature
in European country is smoother than a piecewiselinear function
(Bessec \& Fouquau, 2008), in part due to the nonhomogeneous
@@ 56,48 +35,13 @@ model as simple as possible using the above linear basis.
In addition, the behaviour of the consumers differs
significantly for the week days, Saturdays and Sundays.
Thus, the model is factorized for all three type of days,
resulting in a total of :math:`3 \times 3 = 9` coefficients
+resulting in a total of `9` coefficients
to be adjusted, for each region.
The thresholds :math:`T_H` and :math:`T_C` are two
+The thresholds are two
hyperparameters chosen as constant for all regions and every type of day.
Our strategy is to perform a gridsearch on these thresholds
and to chose their best values via crossvalidation.
Bayesian ridge regression
=========================

For the clarity of the presentation, we focus on a single
region with demand :math:`D` and air temperature :math:`T`
and consider only one type of day.
We use the *scikitlearn* (Buitinck *et al.*, 2013)
opensource implementation of the Bayesian ridge regression of MacKay (1992).
Let :math:`\mathbf{w} = (w_H, w_C, w_0)`
be the vector of basis coefficients in :eq:`demandFunction`.

The noise :math:`\epsilon` in :eq:`demandModel` is assumed to be
Gaussian with zero mean and variance :math:`\alpha`
so that the likelihood
:math:`p(D  T, \mathbf{w}, \alpha)` follows a Gaussian distribution
with mean :math:`X \mathbf{w}` and variance :math:`\alpha`.
The parameter :math:`\alpha` is taken as a random variable.
The prior for :math:`\mathbf{w}` is given by a spherical Gaussian

.. math::
 p(\mathbf{w}  \lambda) = \mathcal{N}(0, \lambda^{1} \mathbf{I}),

where :math:`\lambda` is a random variable and
:math:`\mathbf{I}` is the identity matrix.

The priors over :math:`\alpha` and :math:`\lambda`
are chosen to be gamma distributions,
the conjugate prior for the precision of the Gaussian.
The parameters :math:`\mathbf{w}`, :math:`\alpha` and :math:`\lambda`
are estimated jointly during the fit of the model.
The remaining hyperparameters are the parameters of the gamma priors
over :math:`\alpha` and :math:`\lambda`.
These are chosen to be *noninformative*.
The parameters are estimated by maximizing the *marginal log likelihood*.

.. seealso::
* MacKay, D.J.C., 1992, Bayesian Interpolation, *Neural Computation*, 4, 415447.
diff git a/doc/source/user_guide/generation.rst b/doc/source/user_guide/generation.rst
index 090cf365eb4151eb88f554c8a78edf292e2784c0..739bb65927e53f20767fa297ea16dc42966ab9ff 100644
 a/doc/source/user_guide/generation.rst
+++ b/doc/source/user_guide/generation.rst
@@ 27,26 +27,17 @@ A specific type of wind turbine is selected, the Siemens SWT2.3MW101m.
Hub height correction
.....................
In order the compute the wind magnitude :math:`V(z)` at the hub height :math:`z = 101` m,
the wind magnitude :math:`V(z_0)` at :math:`z_0 = 10` m taken from the CORDEX/HyMex
+In order the compute the wind magnitude at the hub height,
+the wind magnitude taken from the CORDEX/HyMex
simulations is used assuming that the windspeed vertical profile
follows a power law (Justus & Mikhail, 1976), i.e.

.. math::
 V(z) = V(z_0) \left(\frac{z}{z_0}\right)^\alpha,

where value of the exponent :math:`\alpha` is chosen as :math:`1/7`, by default.
+follows a power law (Justus & Mikhail, 1976) where value of the exponent is chosen as `1/7`, by default.
Before to compute the production, two corrections are applied to the wind magnitude.
Air density correction
......................
Deviations of the air density :math:`\rho` from the standard density
:math:`\rho_0` for which the power curve is provided are taken into account
by multiplying the the wind magnitude by a factor

.. math::
 \left(\frac{\rho}{\rho_0}\right)^{1/3}.
+Deviations of the air density from the standard density
+`\rho_0` for which the power curve is provided are taken into account.
This factor is derived by assuming that the wind power is proportional
to the density times the velocity cube. The wind magnitude is corrected
@@ 54,25 +45,16 @@ rather than the wind production, in order to shift the power curve horizontally
rather than scale it vertically, thus respecting the various production regimes.
The air density is computed by applying the ideal gas law to
the atmospheric temperature :math:`T_a`, pressure :math:`P`
and specific humidity :math:`\omega` at the surface taken from the
CORDEX/HyMex data (:py:func:`e4clim.production.getAirDensity`), i.e.

.. math::
 \rho(T_a, P, \omega) = \frac{M_d M_v}{M_v + \omega (M_d  M_v)} \frac{P}{R T_a},

where :math:`M_d` and :math:`M_v` are respectively the molar mass of dry air
and of water vapor, and :math:`R` is the universal gas constant.
+the atmospheric temperature, pressure
+and specific humidity at the surface taken from the
+CORDEX/HyMex data (:py:func:`e4clim.production.getAirDensity`).
Intraday variability correction
................................
In order to take into account that the wind production
is computed from the dailymean windmagnitudes rather than instantaneous values,
the dailymean wind speed is multiplied by a factor

.. math::
 \left(\frac{6}{\pi}\right)^{1/3}.
+the dailymean wind speed is multiplied by a factor.
This coefficient is derived by assuming that the intraday windmagnitude probability
distribution follows a Rayleigh distribution with a scaling factor corresponding to
@@ 107,25 +89,17 @@ the diurnal cycle should be taken into account.
Yet, hourly solar irradiance is usually not available
from the CORDEX/HyMex dataset.
Thus, it is assumed that the clearness index :math:`K_T = I / I_0`,
+Thus, it is assumed that the clearness index,
defined a the ratio between
the surface horizontal irradiance :math:`I`
and the extraterrestrial irradiance :math:`I_0`,
+the surface horizontal irradiance
+and the extraterrestrial irradiance,
is constant throughout the day.
See :py:func:`e4clim.production.getClearnessIndex`.
Then, the hourly extraterrestrial horizontal irradiance is computed
using the Fourier series model of Spencer (1971) to compute
the extraterrestrial radiation :math:`G_0^N` from the calendar data
(:py:func:`e4clim.production.getETRadiation`), i.e.

.. math::
 G_{0}^N = S (1.000110
 + 0.034221 \cos(B) + 0.001280 \sin(B)
 + 0.000719 \cos(2 B) + 0.000077 \sin(2 B)),

where :math:`B = 2 \pi (d  1) / 365` is the phase of the seasonal cycle
and :math:`S = 1367 (W/m2)` is the solar constant.
+the extraterrestrial radiation from the calendar data
+(:py:func:`e4clim.production.getETRadiation`).
See :py:func:`e4clim.production.getETHorizontalIrradiance`.
The hourly extraterrestrial horizontal irradiance is then multiplied
@@ 143,89 +117,49 @@ an approximation of the hourly surface horizontal irradiance.
Computing the tilted horizontal irradiance
..........................................
The tilted horizontal irradiance :math:`I_T` for an anisotropic sky
+The tilted horizontal irradiance for an anisotropic sky
is divided (Duffie & Beckman, 2013)
into a direct (:math:`I_{d, T}`), a diffuse (:math:`I_{r, T}`)
and a reflected component (:math:`I_{g, T}`),

.. math::
 I_T = I_{d, T} + I_{r, T} + I_{g, T}
+into a direct , a diffuse
+and a reflected component.
By default, the array is assumed to face South and to be
tilted by an angle :math:`\beta` equal to the latitude of the array.
+tilted by an angle equal to the latitude of the array.
Arbitrary values and tracking options are, however, also possible.
See :py:func:`e4clim.production.getTiltedIrradiance`.
Direct and diffuse horizontal irradiance
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
The surface horizontal irradiance :math:`I`
is devided into a direct (:math:`I_b`) and a diffuse component (:math:`I_d`).
+The surface horizontal irradiance
+is devided into a direct and a diffuse component.
The ratio between the diffuse and the total horizontal irradiance
is parametrized using the model of Reindl \emph{et al.} (1990a), i.e.

.. math::
 \frac{I_d}{I} =
 \begin{cases}
 1.4  1.749 K_T + 0.177 \sin{e} & \mathrm{if} \quad 0.3 \le K_T \le 0.78 \\
 1.02  0.254 K_T + 0.0123 \sin{e} & \mathrm{if} \quad K_T < 0.3 \\
 0.486 K_T  0.182 \sin{e} & \mathrm{if} \quad K_T > 0.78,
 \end{cases}

where :math:`e` is the elevation angle of the sun.
This ratio is further clipped between 0.1 and 0.97
+is parametrized using the model of Reindl \emph{et al.} (1990a).
+This ratio is further clipped between 0.1 and 0.97.
See :py:func:`e4clim.production.getDiffuseHorizontalIrradiance`.
A threshold :math:`e_\mathrm{cut}` may be defined
in order to account for the obstruction of the
direct horizontal irradiance :math:`I_b = I  I_d`
by objects at solar elevations below :math:`e_\mathrm{cut}`, i.e.

.. math::
 I_b = 0 \quad \mathrm{if} \quad e < e_\mathrm{cut}.

+A threshold may be defined in order to account for the obstruction of the
+direct horizontal irradiance by objects at solar elevations below.
See :py:func:`e4clim.production.getDirectHorizontalIrradiance`.
Direct tilted irradiance
,,,,,,,,,,,,,,,,,,,,,,,,
The transposition factor :math:`R_b` from the horizontal to the tilted
direct irradiance is classically defined as

.. math::
 R_b &= \frac{\cos{\theta}}{\cos{\theta_z}} \\
 \mathrm{so~that~} I_{b, T} &= R_b I_b,

where :math:`\theta_z` is the zenith angle and :math:`\theta` is
the angle of incidence of sun beams on the tilted array.
+The transposition factor from the horizontal to the tilted
+direct irradiance is classically defined.
It is finally ensured that the transposition fact is zero
whenever :math:`I_b` is zero, or when the sun is behind the tilted array.
+whenever the irradiance is zero, or when the sun is behind the tilted array.
See :py:func:`e4clim.production.getDirectTiltedIrradiance`.
Diffuse tilted irradiance
,,,,,,,,,,,,,,,,,,,,,,,,,
The model by Reindl \emph{et al.} (1990b) is used
to compute the diffuse tilted irradiance, i.e.

.. math::
 I_{d, T} = I_d \left\{(1  A_i) \left(\frac{1 + \cos \beta}{2}\right)
 \left[1 + f \sin^3\left(\frac{\beta}{2}\right)\right] + A_i R_b\right\},

where :math:`f = \sqrt{I_b / I}` and :math:`A_i = I_b / I_0` is the anisotropy index,
i.e. the ratio between the direct horizontal irradiance
and the extraterrestrial horizontal irradiance
+to compute the diffuse tilted irradiance.
See :py:func:`e4clim.production.getDiffuseTiltedIrradiance`.
Reflected tilted irradiance
,,,,,,,,,,,,,,,,,,,,,,,,,,,
The reflected part of the irradiance is given by

.. math::
 I_{g, T} = I \gamma_g \left(\frac{1  \cos \beta}{2}\right),

where :math:`\gamma_g` is the surface albedo, with a default value of 0.2.
See :py:func:`e4clim.production.getReflectedTiltedIrradiance`.
.. seealso::
@@ 245,60 +179,18 @@ See :py:func:`e4clim.production.getReflectedTiltedIrradiance`.
Photovoltaic array production
.............................
The electricity production from the photovoltaic array reads
+See :py:func:`e4clim.production.getSolarPower`.
.. math::
 P_\mathrm{PV} = n~A~\eta_\mathrm{e}~\eta_\mathrm{c}~I_T,
+The cell efficiency depends on the cell temperature
+and is computed via the thermal model in Duffie & Beckman (2013).
+If the nominal power and the module area
+are given instead of the reference cell efficiency.
where :math:`n` is the number of modules, :math:`A` their area,
:math:`\eta_\mathrm{e}` is the overall efficiency of the electrical installation,
:math:`\eta_\mathrm{c}` is the efficiency of the cells and :math:`I_T`
is the global tilted irradiance.
See :py:func:`e4clim.production.getSolarPower`.
+The cell temperature is computed from the air temperature,
+the global tilted radiation
+(taken here equal to the hourly irradiance)
+and the wind magnitude.
The cell efficiency :math:`\eta_c` depends on the cell temperature :math:`T_c`
and is computed via the thermal model in Duffie & Beckman (2013), i.e.

.. math::
 \eta_c(T_c) = \eta_{c0} (1  \beta_\mathrm{ref} (T_c  T_{c, \mathrm{ref}})).

where the temperature coefficient :math:`\beta_\mathrm{ref}`
has a default value of :math:`0.004~\mathrm{K}^{1}`,
valid for crystalline silicon modules (Notton \emph{et al.}, 2005).
The reference cell temperature :math:`T_{c, \mathrm{ref}}`
has a value of :math:`25^\circ` C by default.
If the nominal power :math:`P_N` and the module area :math:`A`
are given instead of the reference cell efficiency :math:`\eta_{c0}`,
the latter is computed from

.. math::
 \eta_{c0} = \frac{P_N}{G_{\mathrm{STC}} A},

where :math:`G_{\mathrm{STC}} = 1000~\mathrm{Wm}^{2}`
is the light intensity under Standard Test Conditions (STC).
By default, we set :math:`P_N = 250~\mathrm{Wm}^{2}`
and :math:`A = 1.675`.

The cell temperature is computed from the air temperature :math:`T_a`,
the global tilted radiation :math:`G_T`
(taken here equal to the hourly irradiance :math:`I_T`)
and the wind magnitude :math:`V` as

.. math::
 T_c(T_a, G_T, V) = T_a + \frac{G_T}{G_{T, \mathrm{NOCT}}}
 L(V) (T_{c, \mathrm{NOCT}}  T_{a, \mathrm{NOCT}}),

where the air temperature at NOCT has a value of :math:`20^\circ` C
and the cell temperature at Nominal Operating Cell Temperature (NOCT)
:math:`T_{c, \mathrm{NOCT}}` has a value of :math:`46^\circ` C by default.
The factor :math:`L(V)` depends on the wind speed :math:`V`
and is given by

.. math::
 L(V) = \frac{9.5}{5.7 + 3.8~V}.

If the wind speed is not provided, it is assumed to be constant
and equal to :math:`1~\mathrm{ms^{1}}`.
See :py:func:`e4clim.production.getCellEfficiency`.
.. seealso::
diff git a/doc/source/user_guide/optimization.rst b/doc/source/user_guide/optimization.rst
deleted file mode 100644
index d3bb2ebaceeac43dca64b4eaa4e39365b848cf83..0000000000000000000000000000000000000000
 a/doc/source/user_guide/optimization.rst
+++ /dev/null
@@ 1,231 +0,0 @@
**************************************
Spatial and technological optimization
**************************************

Objective
=========

The objective is to optimise the penetration rate of renewable energies
in the Italian electricity mix by using a spatial and technological diversification
between photovoltaic and wind power across different grid regions.

Meanvariance analysis
======================

The geographical and technological diversification of the power plants
is carried out in accordance with an analysis that is directly derived
from Markowitz's modern portfolio theory.

A meanvariance analysis is the process of assessing the risk (variance)
with respect to a expected result (average).
By examining the penetration rate (mean) and variance of a
space and technology assembly of renewable power plants, it is possible to make
more efficient installation choices  by seeking the lowest variance for a desired rate or by
looking for the highest rate for a given variance.

Each combination of renewable installations can be represented by a meanvariance graph.
For each mean, there is a combination that minimizes variance.
Conversely, for At each level of variance, a portfolio can be found
that maximizes the expected average.
All of these portfolios are called efficient frontier or Markowitz frontier.
Border points are by definition suboptimal,
and will not be of interest to a rational investor.
The area above the border cannot be reached.

Method
======

A python minimization function, provided by the Scipy package, is used to
determine this efficient boundary.
For a given mean interval (between 0 and 50\%),
the set of combination of renewable facilities minimizing variance
is determined.
Using this minimization function, we implemented three methods
of increasing complexity to perform the meanvariance analysis:

* The first procedure is the most ideal.
 The electricity generated at a point is assumed to be immediately available to meet the
 overall Italian demand, regardless of any distance and transport constraints (Eq. 1 \& 2).
* The second includes a notion of electricity network.
 Local production in a region must first be used to provide local demand.
 If production is stronger than demand, electricity may be exported
 to the extent that the capacity of the power lines permits.
 Otherwise the electricity is lost (Eq. 3).
* Finally, we add a financial component on the profitability of the facilities
 renewable (E. 4).

Global optimization


In the first case, the following objective function is minimized,

.. math::
 \sum_{i = 1}^N \sum_{j = 1}^N \omega_i \omega_j \sigma_{ij}
 + 1e6 \left(\left( \sum_{i = 1}^N \omega_i \mu_i\right)  \mu_{\mathrm{Target}}\right),

where :math:`\mu_i(t)` is the average penetration level per region and technology
for an installed MW compared to the national electricity demand,
:math:`N` is the number of regions by the number of available technologies
(2, photovoltaic and wind), :math:`\omega` is the installed capacity in MW,
:math:`\mu_\mathrm{Target}` is the target penetration rate
and :math:`\sigma` is the covariance of the renewable penetration rate.
The average penetration level per region and technology is given by,

.. math::
 \mu_i = \frac{1}{T} \sum_{t = 1}^T \frac{P_i(t)}{\sum_{i = 1}^N D_i(t)},

where :math:`P_i(t)` is the daily production time series
by region and technology :math:`i`
and :math:`D_i(t)` is the power consumption.

.. seealso::
 solve_frontier_portfolio

Regional optimization


In the second optimization problem, the interconnection between regions
is taken into account. The following constraints are thus added,

.. math::
 \mu_i = \frac{1}{T} \sum_{t = 1}^T \frac{P_i(t)}{D_i(t)},

 P_i(t) \le D_i(t) + T_i

where :math:`T_i` is the electrical transmission capacity of a region.

.. seealso::
 solve_frontier_portfolio_constraint

Return on expenditures optimization


The last optimization problem adds the following constraint
ensuring a return on expenditures,

.. math::
 \sum_{i = 1}^N \sum_{j = 1}^N \omega_i \omega_j \sigma_{ij}
 + 1e6 \left(\left( \sum_{i = 1}^N \omega_i \mu_i\right)  \mu_{\mathrm{Target}}\right)
 + \frac{1}{N} \sum_{i = 1}^N
 \frac{\omega_i (\mathrm{CAPEX}_i + 24 \mathrm{OPEX}_i)}{
 \frac{\omega_i p_\mathrm{Elec}}{T} \sum_{t = 1}^T \mu_i(t) D_i(t)},

where :math:`\mathrm{CAPEX}_i` is the initial capital expenditure
whereas :math:`\mathrm{OPEX}_i` is the operating expenditure
(current maintenance and operating costs),
:math:`p_\mathrm{Elec}` is the electricity price.

.. seealso::
 solve_frontier_portfolio_roe

Files and programs
==================

* `Opt_portfolio.py`: Calculates the efficient frontier for the 3 optimization levels. Store them
 results for configurations with minimum variance and maximum mean.
* `Ratio_SpecialDays_Script.py`: Calculates the percentage of situation at risk for the 3 levels
 optimizations and following distributions ranging from photovoltaic to wind energy.
 The situations at risk are when,

 1. the production does not reach the level of national demand,
 2. the penetration rate of renewables exceeds 40% of the electricity mix.

* `DrawDistrib.py`: Draw the optimization results on a map.


Alexis
======

Let :math:`N` denote the number of regions in the domain.
Two technologies are considered, photovoltaic solar panels
and wind turbines.
Let :math:`\mathbf{k} = (i, j), 1 \le i \le N, 1 \le j \le 2`,
be a multiindex with :math:`i` being the regional index
and :math:`j` being the technological index.

The capacity factors for each region and technology are defined by

.. math::
 \eta_\mathbf{k}(t) = \frac{P_\mathbf{k}(t)}{P^N_\mathbf{k}}

Let :math:`D_T(t)` denote the total demand at time :math:`t`, i.e.

.. math::
 D_T(t) = \sum_i D_i(t),

where :math:`D_i(t)` is the demand of region :math:`i`.
Let :math:`(\omega_\mathbf{k})_\mathbf{k}` be the vector
of RES capacity (in MW) per region and technology and

.. math::
 \mu_\mathbf{k}(t) = \frac{\eta_\mathbf{k}(t)}{D_T(t)},

be the penetration rates per region and technology at time :math:`t`, i.e.
the fraction of the total demand satisfied by the RES production in a given region.
The total penetration of RES in the electricity mix :math:`\mu_T(t)`
is then given by

.. math::
 \mu_T(t) = \sum_\mathbf{k} \omega_\mathbf{k} \mu_\mathbf{k}(t).


Our objective is then to maximize the mean total penetration

.. math::
 \mathbb{E}[\mu_T] = \sum_\mathbf{k} \omega_\mathbf{k} \mathbb{E}[\mu_\mathbf{k}],

while minimizing some measure of the risk.
The risk is usually defined as the variance of :math:`\mu_T(t)`, i.e. as

.. math::
 \mathbb{V}[\mu_T] &= \sum_{\mathbf{k}, \mathbf{l}} \omega_\mathbf{k} \sigma_{\mathbf{k} \mathbf{l}} \omega_\mathbf{l},

where the matrix :math:`\Sigma` is such that

.. math::
 \sigma_{\mathbf{k} \mathbf{l}} = \mathbb{E}[(\mu_\mathbf{k} \mu_\mathbf{l}  \mathbb{E}[\mu_\mathbf{k} \mu_\mathbf{l}])^2]

However, in order to value the covariance of the regional RES production with respect to the regional demand,
rather than with respect to the aggregated demand, we prefer
define the risk by replacing the matrix :math:`\Sigma`
with the matrix :math:`Q` defined such that

.. math::
 q_{\mathbf{k} \mathbf{l}} = \mathbb{E}[(\nu_\mathbf{k} \nu_\mathbf{l}  \mathbb{E}[\nu_\mathbf{k} \nu_\mathbf{l}])^2],

with the vector :math:`(\eta_\mathbf{k})_\mathbf{k}`
giving the fraction of the regional demand satisfied by the RES production, i.e.

.. math::
 \nu_\mathbf{k}(t) = \frac{\eta_\mathbf{k}(t)}{D_i(t)}.

Thus, we want to minimize the risk

.. math::
 \sum_{\mathbf{k}, \mathbf{l}} \omega_\mathbf{k} q_{\mathbf{k} \mathbf{l}} \omega_\mathbf{l}.

Note, however, that this risk no longer coincides with the variance of the total penetration rate.

To solve this biobjective problem, we minimize the risk while constraining
the mean penetration rate to some target and repeat this operation for different targets.
That is we solve the problem

.. math::
 \min_{\omega} &\sum_{\mathbf{k}, \mathbf{l}} \omega_\mathbf{k} q_{\mathbf{k} \mathbf{l}} \omega_\mathbf{l} \\
 \mathrm{subject~to} \quad &\sum_\mathbf{k} \omega_\mathbf{k} \mathbb{E}[\mu_\mathbf{k}] = U.

We also impose the following upper and lower bounds on the capacities

.. math::
 0 \le \omega_\mathbf{k} \le \overline{\omega},

where :math:`\overline{\omega}` is a maximum capacity.

To avoid congestion situations within regions, we also add the following constraint

.. math::
 \omega_\mathbf{k} \eta_\mathbf{k}(t) \le D_i(t) + T_i,

where :math:`T_i` is the transmission capacity from region :math:`i`
to any other region.

diff git a/e4clim/actuator/pv.py b/e4clim/actuator/pv.py
index 7413323b1b37900e2cf18cc41b0fbc5a5d999f47..d33d625bee284cfa4887bf516b6c5e02f3bc2f59 100644
 a/e4clim/actuator/pv.py
+++ b/e4clim/actuator/pv.py
@@ 287,7 +287,7 @@ def get_pv_power(global_tilted_surf, temp_a=None, wind_speed=1., area=1.675,
temp_ref. If None, compute it from the nominal power. Default is None.
:param thermal_loss: Efficiency loss per degree Kelvin
away from the reference temperature. Default is
 :math:`0.0041~\mathrm{K}^{1}`.
+ `0.0041 K1`.
:param temp_ref: Reference temperature (K). Default 298.15 K.
:param temp_cell_noct: Cell temperature at the
Nominal Operating Cell Temperature (NOCT) (K).
@@ 361,7 +361,7 @@ def get_cell_efficiency(
temp_ref. If None, compute it from the nominal power. Default is None.
:param thermal_loss: Efficiency loss per degree Kelvin
away from the reference temperature. Default is
 :math:`0.0041~\mathrm{K}^{1}`.
+ `0.0041 K1`.
:param temp_ref: Reference temperature (K). Default 298.15 K.
:param temp_cell_noct: Cell temperature at _noct (K). Default is 319.15 K.
:param nominal_power: Nominal power of the module (W). See note.
@@ 410,11 +410,11 @@ def get_cell_efficiency(
temp_cell = temp_ref
else:
# Compute temperature dependent efficiency of module
 # Thermal loss ration: :math:`U_{L, _noct} / U_L`
+ # Thermal loss ration: U_{L, _noct} / U_L
loss_ratio = 9.5 / (5.7 + 3.8 * wind_speed)
# Cell temperature (K). See note.
 # :math:` au \alpha / U_L = (T_{c, _noct}  T_{a, _noct}) / G_{T, NOCT}`
+ # au \alpha / U_L = (T_{c, _noct}  T_{a, _noct}) / G_{T, NOCT}
temp_cell = (temp_a + (global_tilted_surf / global_tilted_surf_noct)
* loss_ratio * (temp_cell_noct  temp_a_noct))
diff git a/e4clim/actuator/wind.py b/e4clim/actuator/wind.py
index 846f18e025be18a6f1803f3d87ebaa2a830b1747..b4ed70b5c3611b769ada6e6c8cc2cadbcdca98f4 100644
 a/e4clim/actuator/wind.py
+++ b/e4clim/actuator/wind.py
@@ 287,9 +287,7 @@ class Actuator(ExtractorBase):
:param rho0: Reference air density (kg/m3) for which the
power curve was obtained. Default is 1.225 (kg/m3).
:param mean_speed2mean_gen: If True, :py:obj:`speed`
 is taken as the mean of a Rayleigh distribution
 so that a factor :math:`(6 / \pi)^{1/3}` is applied to the
 mean wind speed to get the corresponding mean wind power.
+ is taken as the mean of a Rayleigh distribution.
Default is False.
:type speed: (sequence of) :py:class:`float`
:type rho: (sequence of) :py:class:`float`
diff git a/e4clim/optimizer/res_mean_risk.py b/e4clim/optimizer/res_mean_risk.py
index a1a317f95e1dfba09101aff957b4c4facad341c6..20b109d0dfe9ed5eeb67c2acb0eadea84d05862f 100644
 a/e4clim/optimizer/res_mean_risk.py
+++ b/e4clim/optimizer/res_mean_risk.py
@@ 1734,8 +1734,8 @@ def fun_regional_transmission(x, *args):
omega_t_wind = omega_t[:, n_reg:]
cf_t_pv = args[0][:, :n_reg]
cf_t_wind = args[0][:, n_reg:]
 # :math:`T_i+D_i(t)\omega_i \eta_i(t)\omega_{i + n_reg}
 # \eta_{i + n_reg}(t)`
+ # T_i+D_i(t)\omega_i \eta_i(t)\omega_{i + n_reg}
+ # \eta_{i + n_reg}(t)
y = (trans_t + dem_t  omega_t_pv * cf_t_pv 
omega_t_wind * cf_t_wind).flatten()
@@ 1962,10 +1962,10 @@ def get_bounds_cstr(lb=None, ub=None, dim=None, sparse=False):
:type dim: int
:type sparse: bool
 :returns: A tuple containing matrix :math:`G` and vector
 :math:`h` defining constraint :math:`G x \le h`.
+ :returns: A tuple containing matrix `G` and vector
+ `h` defining constraint `G x <= h`.
:rtype: tuple of :py:class:`numpy.array`,
 or :py:obj:`cvxspmatrix` and :py:obj:`cvxmatrix`.
+ or :py:class:`cvxopt.spmatrix` and :py:class:`cvxopt.cvxmatrix`.
.. note:: Constraint is adapted to the quadprod module.
To use the CVXOPT module, the oposite should be taken.
@@ 2000,7 +2000,7 @@ def get_bounds_cstr(lb=None, ub=None, dim=None, sparse=False):
G.append(np.eye(dim))
h.append(uba)
 # Matrix G and vector h for :math:`G x \le h`
+ # Matrix G and vector h for `G x <= h`
h = np.concatenate(h)
if sparse:
v = np.concatenate(v)
@@ 2023,9 +2023,9 @@ def get_weighted_sum_cstr(weights, target_sum):
:type weights: :py:class:`numpy.array`
:type target_sum: float
 :returns: A tuple containing matrix :math:`A` and vector
 :math:`b` defining constraint :math:`A x = b`.
 :rtype: tuple of :py:class:`numpy.array`.
+ :returns: A tuple containing matrix `A` and vector
+ `b` defining constraint `A x = b`.
+ :rtype: tuple of :py:class:`numpy.ndarray`.
"""
return (np.expand_dims(weights, 0), np.array([target_sum]))
@@ 2103,17 +2103,12 @@ def distribute_cv_square(power_mismatch, r):
all nodes but the first one to be positive.
Without this bound, one would have:
 * For :math:`r = 0`, the conventional production is uniformly
+ * For `r = 0`, the conventional production is uniformly
distributed.
 * For :math:`r = 1`, the conventional production compensates
+ * For `r = 1`, the conventional production compensates
the mismatch between the RES produciton and the demand,
i.e. there is no transmission.
 * For :math:`0 \le r \le 1`, the conventional production is
 given by

 .. math::
 (1  r) \sum_{i = 1} ^ N \frac{P ^ {\mathrm{_res}}_i  dem_i}{N}
  r(P ^ {\mathrm{_res}}_i  dem_i)
+ * For `0 <= r \<= 1`, the conventional production by a mix of the two.
* Since may(for each time step) small
(the size given by the number of regions)
@@ 2267,13 +2262,7 @@ def distribute_cv_mean_square(power_mismatch, r, window_length=None):
def quad_prog_active_set(G, c, A, b=None):
 """Bruteforce active set method to solve convex quadratic program

 .. math::

 \min_x x^t G x + c^t x

 \mathrm{subject~to~} A^t x >= b.
+ """Bruteforce active set method to solve convex quadratic program.
:param G: Positive definite matrix.
:param c: Linear vector.
@@ 2289,7 +2278,7 @@ def quad_prog_active_set(G, c, A, b=None):
* Optimization problem solution.
* Lagrange multipliers.
 :rtype: tuple of :py:class:`~numpy.array`
+ :rtype: tuple of :py:class:`numpy.ndarray`
.. note:: This function is used for pedagogical purposes only.
Favor using more advanced solvers (e.g. quadprog or from scipy).
diff git a/e4clim/utils/solar_calendar.py b/e4clim/utils/solar_calendar.py
index f4143babb5a056ec315c75eb2e7a9ee5f02a1760..cfd27abf1185986769332b1661095c071d07e041 100644
 a/e4clim/utils/solar_calendar.py
+++ b/e4clim/utils/solar_calendar.py
@@ 55,7 +55,7 @@ def _compute_if_none(fun):
def _clip_angle(angle):
 """Clip angle so that: math: `\pi ^\\circ <= \\theta < \pi ^\\circ`.
+ """Clip angle so that `pi <= theta < pi`.
:param angle: Angle in radians.
:type angle: (sequence of) :py:class:`float`
@@ 90,7 +90,7 @@ class SolarCalendarComputer(object):
* The *declination* is the angular position of the sun at solar noon
(i.e., when the sun is on the local meridian)
with respect to the plane of the equator, north positive;
 :math:`−23.45^\circ \le \delta \le 23.45^\circ`.
+ `−23.45° <= \delta <= 23.45°`.
* The *hour angle* is the angular displacement of the sun east or west
of the local meridian due to rotation of the earth on its axis at
15 degrees per hour; morning negative, afternoon positive.
@@ 125,7 +125,7 @@ class SolarCalendarComputer(object):
.. warning::
* All attribute variables (e.g. :py:attr:`hour_angle`)
 are given in radians, between :py:math:`\pi` and :py:math:`\pi`,
+ are given in radians, between `pi` and `pi`,
and all computations are performed
in radians. To get angles in the same units as those
provided to constructor, call :py:meth:`get_angle`.
@@ 197,7 +197,7 @@ class SolarCalendarComputer(object):
@property
@_compute_if_none
def mean_longitude(self):
 """Get mean longitude :py:math:`\pi <= mean_lon <= \pi`."""
+ """Get mean longitude `pi <= mean_lon <= pi`."""
return 4.89495 + 0.01720279 * self.ndays_since
@property