opti prod

opti_prod.py

@@ -0,0 +1,60 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import PyQt5 as qt
+
+ptot = 1000 # [MW]
+
+def C1(x):
+    return 30*x + 0.01*x**2
+
+def C2(x):
+    return 20*x + 0.02*x**2
+
+def f(x):
+    return C1(x[0]) + C2(x[1]) + x[-1] * (ptot - x[0] - x[1])
+
+def grad(f, x, h=1e-4):
+    res = []
+    for i in range(len(x)):
+        delta = f(x[:i] + [x[i] + h / 2] + x[i+1:]) - f(x[:i] + [x[i] - h / 2] + x[i+1:])
+        res += [delta / h]
+    return res
+
+def norm(x):
+    n = 0
+    for d in x:
+        n += d**2
+    return np.sqrt(n)
+
+def g(x):
+    return norm(grad(f, x))
+
+def minize(f, x0, h=1e-4, step=1e-1, tol=1e-8, N=1e4, echo=False):
+    x = x0
+    g = grad(f, x, h)
+    print(g)
+    n = 0
+    prev = norm(g) + 2*tol
+    print(prev, norm(g), abs(norm(g) - prev))
+    while abs(norm(g) - prev) > tol:
+        n += 1
+        prev = norm(g)
+        for i in range(len(x)):
+            x[i] -= g[i] * step
+        g = grad(f, x, h)
+        if (n % 100 == 0) and echo:
+            print("Itération ", n)
+            print("norm(g) = ", norm(g))
+            print("prev = ", prev)
+            print("x = ", x)
+            print("g = ", g)
+        if n > N:
+            return x
+    
+    return x
+
+#print(f([500, 500, 40]))
+#print(f([450, 450, 35]))
+#print(C1(500))
+#print(C2(500))
+print(minize(g, [0, 0, 0]))