Antares Simulator
Power System Simulator
math.hxx
1 /*
2  * Copyright 2007-2025, RTE (https://www.rte-france.com)
3  * See AUTHORS.txt
4  * SPDX-License-Identifier: MPL-2.0
5  * This file is part of Antares-Simulator,
6  * Adequacy and Performance assessment for interconnected energy networks.
7  *
8  * Antares_Simulator is free software: you can redistribute it and/or modify
9  * it under the terms of the Mozilla Public Licence 2.0 as published by
10  * the Mozilla Foundation, either version 2 of the License, or
11  * (at your option) any later version.
12  *
13  * Antares_Simulator is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * Mozilla Public Licence 2.0 for more details.
17  *
18  * You should have received a copy of the Mozilla Public Licence 2.0
19  * along with Antares_Simulator. If not, see <https://opensource.org/license/mpl-2-0/>.
20  */
21 
22 #include <cmath>
23 
24 #include <yuni/yuni.h>
25 
26 #include <antares/logs/logs.h>
27 #include <antares/study/study.h>
28 #include "antares/solver/ts-generator/xcast/xcast.h"
29 
30 #include "constants.h"
31 
32 using namespace Yuni;
33 
34 namespace Antares::TSGenerator::XCast
35 {
36 /*
37 ** This function is no longer used but remains for code clarity
38 **
39 float Minimum(float a,float b,float g,float d,int l)
40 {
41  switch(l)
42  {
43  case 1: return 0;
44  case 2: return 0;
45  case 3: return 0;
46  case 4: return 0;
47  case 5: return 0;
48  }
49  return 0.f;
50 }
51 */
52 
53 static float maximum(float a, float b, float g, float d, int l)
54 {
55  switch (Data::XCast::Distribution(l))
56  {
57  case Data::XCast::dtBeta:
58  case Data::XCast::dtUniform:
59  case Data::XCast::dtNormal:
60  {
61  return d - g;
62  }
63  case Data::XCast::dtWeibullShapeA:
64  {
65  // on regle la borne pour que le residu sur la fonction de repartition :
66  // y= 1-exp(-(x/b)^a) soit < NEGLI_WEIB
67  return (b * (float(pow(log(1. / NEGLI_WEIB), double(1. / a))))) - g;
68  }
69  case Data::XCast::dtGammaShapeA:
70  {
71  // on regle la borne pour que la fonction de distribution :
72  // y = b*((bx)^a-1)*exp(-bx)/gamma_euler(a) soit < NEGLI_GAMM
73  // on regle la borne ou une valeur pessimiste en considerant que la
74  // fonction ne diminue pas plus rapidement qu'une exponentielle
75  // (alors qu'on impose a >= 1)
76  // Since 3.5, B' = 1/B
77  return ((float)(log(1. / NEGLI_GAMM) * a * b)) - g;
78  }
79 
80  case Data::XCast::dtMax:
81  case Data::XCast::dtNone:
82  return 0.f;
83  }
84  return 0.f;
85 }
86 
92 static float esperance(float a, float b, float g, float d, int l)
93 {
94  // note : valeur exprimee en valeur relative par rapport au minimum
95  switch (Data::XCast::Distribution(l))
96  {
97  case Data::XCast::dtBeta:
98  return (a * (d - g)) / (a + b);
99 
100  case Data::XCast::dtUniform:
101  return (d - g) / 2.f;
102 
103  case Data::XCast::dtNormal:
104  return 6.f * b;
105 
106  case Data::XCast::dtWeibullShapeA:
107  return b * (float)XCast::GammaEuler(1. + 1. / double(a));
108 
109  case Data::XCast::dtGammaShapeA:
110  return a * b;
111 
112  case Data::XCast::dtMax:
113  case Data::XCast::dtNone:
114  return 0.f;
115  }
116  return 0.f;
117 }
118 
120 static float standard(float a, float b, float g, float d, int l)
121 {
122  switch (Data::XCast::Distribution(l))
123  {
124  case Data::XCast::dtBeta:
125  {
126  float s = a + b;
127  s = s * s * (s + 1.f);
128  s = a * b / s;
129  s = (float)sqrt(s);
130  s = (d - g) * s;
131  return s;
132  }
133  case Data::XCast::dtUniform:
134  {
135  return (d - g) / (3.464101615f /* float(sqrt(12.)) */);
136  }
137  case Data::XCast::dtNormal:
138  {
139  return b;
140  }
141  case Data::XCast::dtWeibullShapeA:
142  {
143  const float x = (float)(XCast::GammaEuler(1. + 2. / double(a))
144  - pow(XCast::GammaEuler(1. + 1. / double(a)), 2.));
145  if (x < 0.f)
146  {
147  logs.warning()
148  << "TS Generator: error when computing the standard deviation (Weibul Shape A)";
149  return 0.f;
150  }
151  return b * sqrt(x);
152  }
153  case Data::XCast::dtGammaShapeA:
154  {
155  return (float)(sqrt(a) * b);
156  }
157 
158  case Data::XCast::dtMax:
159  case Data::XCast::dtNone:
160  return 0.f;
161  }
162  return 0.f;
163 }
164 
166 static float diffusion(float a, float b, float g, float d, int l, float t, float x)
167 {
168  switch (Data::XCast::Distribution(l))
169  {
170  case Data::XCast::dtBeta:
171  {
172  return (float)sqrt(2. * t * x * ((d - g) - x) / (a + b));
173  }
174  case Data::XCast::dtUniform:
175  {
176  return (float)sqrt(t * x * ((d - g) - x));
177  }
178  case Data::XCast::dtNormal:
179  {
180  return b * (float)sqrt(2. * t);
181  }
182  case Data::XCast::dtWeibullShapeA:
183  {
184  // Warning:
185  // On suppose ici que d a ete correctement initialise avec d=gamma_euler(1+1/a)
186  //
187  assert(b != 0 && "division by zero");
188  const double v = pow(x / b, a);
189  const double w = exp(v);
190 
191  double y = (w - 1.) * double(d);
192  y -= w * XCast::GammaInc(double(1. + 1. / a), v);
193  y /= v;
194  y *= x * b * 2. * t / a;
195  if (y < 0.)
196  {
197  logs.info() << " square diff: " << y << " (Weibul Shape A)";
198  y = 1e-6;
199  }
200  return (float)sqrt(y);
201  }
202  case Data::XCast::dtGammaShapeA:
203  {
204  return (float)sqrt(2. * t * x * b);
205  }
206 
207  case Data::XCast::dtMax:
208  case Data::XCast::dtNone:
209  return 0.f;
210  }
211  return 0.f;
212 }
213 
214 /*
215 ** \brief Verification de l'acceptabilite d'une factorisation LtL
216 **
217 **
218 ** Ce code verifie si LtL est proche de A, A symetrique reelle et L triangulaire
219 ** inferieure qui n'est pas necessairement definie positive.
220 **
221 ** \param L estimation a priori de la racine carrée autoconjuguée de A
222 ** \param A matrice symétrique NxN
223 ** \param N Nombre de processus
224 **
225 ** \return True si la factorisation est correcte a la tolerance pres
226 */
227 /*static bool proma(float** L, float** A, uint N)
228 {
229  float estim_aij;
230  uint i, j, k;
231 
232  for (i = 0; i < N; i++)
233  {
234  for (j = 0; j <= i; ++j)
235  {
236  estim_aij = 0.f;
237  for (k = 0; k < N; ++k)
238  estim_aij += L[i][k] * L[j][k];
239 
240  if (std::abs(estim_aij - A[i][j]) > EPSIMAX)
241  return false;
242  }
243  }
244  return true;
245 }*/
246 
247 static bool verification(float a, float b, float g, float d, int l, float t)
248 {
249  switch (Data::XCast::Distribution(l))
250  {
251  case Data::XCast::dtBeta:
252  {
253  return !(d < g || g < INFININ / 2.f || d > INFINIP / 2.f || t < 0.f || a < 0.f || b < 0.f);
254  }
255  case Data::XCast::dtUniform:
256  {
257  return !(d < g || g < INFININ / 2.f || d > INFINIP / 2.f || t < 0.f);
258  }
259  case Data::XCast::dtNormal:
260  {
261  return !(d < g || g < INFININ / 2.f || d > INFINIP / 2.f || a < INFININ / 2.f
262  || a > INFINIP / 2.f || t < 0.f || b < 0.f);
263  }
264  case Data::XCast::dtWeibullShapeA:
265  {
266  return !(g < INFININ / 2.f || a < 1.f || a > 50.f || t < 0.f || b <= 0.f);
267  }
268  case Data::XCast::dtGammaShapeA:
269  {
270  return !(g < INFININ / 2.f || a < 1.f || a > 50.f || t < 0.f || b <= 0.f);
271  }
272 
273  case Data::XCast::dtMax:
274  case Data::XCast::dtNone:
275  return false;
276  }
277  return false;
278 }
279 
281 static float maxiDiffusion(float a, float b, float g, float d, int l, float t)
282 {
283  switch (Data::XCast::Distribution(l))
284  {
285  case Data::XCast::dtBeta:
286  {
287  return (float)((d - g) * sqrt(t / (2. * (a + b))));
288  }
289  case Data::XCast::dtUniform:
290  {
291  return (float)((d - g) * sqrt(t) / 2.);
292  }
293  case Data::XCast::dtNormal:
294  {
295  return (float)(b * sqrt(2. * t));
296  }
297  case Data::XCast::dtWeibullShapeA:
298  {
299  // pour une loi de Weibull, on prend comme reference le max des diffusions obtenues pour :
300  // x= e= gamma_euler(1+1/a)
301  // x = Maximum
302  // quand a>1, x = b*(((a-1)/a)^(1/a))
303  // la valeur de la diffusion
304 
305  // necessaire si d n'a pas ete initialise au prealable
306  d = float(XCast::GammaEuler(1. + 1. / double(a)));
307 
308  float x = b * d;
309  float m = diffusion(a, b, g, d, l, t, x);
310  float y = diffusion(a, b, g, d, l, t, maximum(a, b, g, d, l));
311  if (y > m)
312  {
313  m = y;
314  }
315  if (a > 1.f)
316  {
317  x = (a - 1.f) / a;
318  x = (float)pow((double)x, (double)(1. / a));
319  x *= b;
320  y = diffusion(a, b, g, d, l, t, x);
321  if (y > m)
322  {
323  m = y;
324  }
325  }
326  return m;
327  }
328  case Data::XCast::dtGammaShapeA:
329  {
330  return diffusion(a, b, g, d, l, t, maximum(a, b, g, d, l));
331  }
332 
333  case Data::XCast::dtMax:
334  case Data::XCast::dtNone:
335  return 0.f;
336  }
337  return 0.f;
338 }
339 
340 } // namespace Antares::TSGenerator::XCast
Antares::Data::XCast::Distribution
Distribution
All available probability distribution.
Definition: xcast.h:61