Merge branch 'develop' of gitlab.mech.kuleuven.be:meco-software/cpp_s…

…plines into develop

examples/global_optimization.py

@@ -0,0 +1,107 @@
+import numpy as np
+import time
+import matplotlib.pyplot as plt
+
+import meco_binaries
+meco_binaries(cpp_splines="develop")
+from splines import *
+
+coef_f_poly = [10,1,4,0,-1]
+f = Polynomial(coef_f_poly)
+
+
+r = np.linspace(-2,2,201)
+f_eval = f.list_eval(r)
+
+plt.subplot(2,1,1)
+plt.plot(r, f_eval)
+
+b = BSplineBasis([-2,2],1, 3)
+
+g = f.transform_to(b)
+
+plt.plot(g.basis().greville(), g.coeff_tensor().flatten())
+opti = OptiSpline()
+
+alpha = opti.variable()
+# opti.subject_to(g<=alpha)  Werkt niet
+opti.subject_to(alpha - g >= 0)
+opti.solver('ipopt')
+opti.minimize(alpha)
+sol = opti.solve()
+
+plt.plot([-2, 2], [sol.value(alpha)] *2)
+
+plt.subplot(2,1,2)
+plt.plot(r, f_eval)
+
+b = BSplineBasis([-2,2],1, 6)
+
+g = f.transform_to(b)
+
+plt.plot(g.basis().greville(), g.coeff_tensor().flatten())
+
+opti = OptiSpline()
+
+alpha = opti.variable()
+# opti.subject_to(g<=alpha)  Werkt niet
+opti.subject_to(alpha - g >= 0)
+opti.solver('ipopt')
+opti.minimize(alpha)
+sol = opti.solve()
+
+plt.plot([-2, 2], [sol.value(alpha)] *2)
+plt.show()
+
+#########################################################
+### simple multi variate global optinmization problem ###
+#########################################################
+grens = 1.25
+knots = [0.5]
+b = BSplineBasis([0, 0, 0.8, grens, grens], 1)
+# b = BSplineBasis([0, grens], 1, 6)
+
+x = Polynomial([0,1],'x')
+y = Polynomial([0,1],'y')
+
+F = x**4 * y**2 +  x**2 * y**4 - 3 * x**2 * y**2 + 1
+
+F = F.transform_to(TensorBasis([ b,b ]))
+
+args = ['x', 'y']
+
+from mpl_toolkits.mplot3d import Axes3D
+
+import matplotlib.pyplot as plt
+from matplotlib import cm
+from matplotlib.ticker import LinearLocator, FormatStrFormatter
+import numpy as np
+
+fig = plt.figure()
+ax = fig.gca(projection='3d')
+
+# Make data.
+X_grid = np.linspace(0,grens, 100)
+Y_grid = np.linspace(0,grens, 100)
+grid = np.meshgrid(X_grid, Y_grid)
+
+Z = F.grid_eval([X_grid, Y_grid])[:,:,0,0]
+
+# Plot the surface.
+X_G, Y_G = np.meshgrid(X_grid, Y_grid)
+surf = ax.plot_surface(X_G, Y_G, Z,
+                       linewidth=0, antialiased=False, alpha=0.5)
+
+X_wire, Y_wire = np.meshgrid(F.basis(0).greville() , F.basis(1).greville())
+print(X_wire)
+print(Y_wire)
+C = F.coeff_tensor()[:,:,0,0]
+ax.plot_wireframe(X_wire, Y_wire, C, linewidth=0.8 , color="black")
+X_G, Y_G = np.meshgrid([0, grens], [0, grens])
+surf = ax.plot_surface(X_G, Y_G, np.ndarray.min(C),
+                       linewidth=0, antialiased=False, alpha=0.3, color="grey")
+# ax.plot_wireframe(F.basis(0).greville() , F.basis(1).greville(), F.coeff_tensor()[:,:,0,0], rstride=10, cstride=10)
+ax.set_xlim(0, grens)
+ax.set_ylim(0, grens)
+ax.set_zlim(-0.4, 1.4)
+plt.show()

examples/simple_controle_example.py

@@ -0,0 +1,84 @@
+import meco_binaries
+meco_binaries(cpp_splines="")
+from splines import *
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib import cm
+from matplotlib.ticker import LinearLocator, FormatStrFormatter
+
+opti = OptiSpline()
+
+class Circle:
+    def __init__(self, center, r):
+        self.x = center[0]
+        self.y = center[1]
+        self.r = r
+
+    def distance_to(self, obj):
+        return (obj[0] - self.x)**2 + (obj[1] - self.y)**2 - self.r **2
+
+    def plot(self, plt):
+        a_ = np.linspace(0,np.pi*2,200)
+        plt.plot(self.r * np.cos(a_) + self.x,self.r * np.sin(a_) + self.y)
+
+obsticals = [Circle((0.3, 0.1), 0.2), Circle((0.5, 0.8), 0.3)]
+
+v_max = 1  #[m/s]
+
+# Define spline
+deg = 3             #degree
+n = 15              #number of knots
+
+basis = BSplineBasis([0, 1], deg, n)
+
+x = opti.Function(basis, [ 2, 1 ])
+v = x.derivative(1)  # time derivative
+
+total_time = opti.variable() #motion time, function of d
+
+# start condictions
+opti.subject_to(x(0) == [0,0])
+
+# end conditions
+opti.subject_to(x(1)  == [1,1.4])
+
+# global constraints
+for o in obsticals:
+    opti.subject_to(o.distance_to(x) >= 0)
+
+opti.subject_to(v[0]**2 + v[1]**2 <= total_time*v_max)
+
+opti.subject_to(total_time >= 0)
+
+# warm starting solver
+opti.set_initial(x, [1, 1.4] * Parameter())
+
+opti.minimize(total_time)
+opti.solver("ipopt",{"ipopt":{"tol":1e-5}})
+
+sol = opti.solve()
+x_sol = sol.value(x)
+v_sol = sol.value(v)
+total_time_sol = sol.value(total_time)
+
+print("total_time_sol : " + str( total_time_sol))
+
+# Make data.
+T_ = np.linspace(0,1,100)
+time = total_time_sol *T_
+x_pos  = x_sol.list_eval(T_)
+fig = plt.figure()
+ax = fig.add_subplot(1, 2, 1)
+ax.plot(x_pos[ :,0 ], x_pos[ :,1 ])
+
+a_ = np.linspace(0,np.pi*2,200)
+
+for o in obsticals:
+    o.plot(ax)
+ax.axis([-0.1, 1.1, -0.1, 1.5])
+ax = fig.add_subplot(1, 2, 2)
+v_abs = (v_sol[0]**2 + v_sol[1]**2) * (1 / total_time_sol)
+ax.plot(total_time_sol * T_, v_abs.list_eval(T_))
+
+ax.axis([0, total_time_sol, 0, 1.2])
+plt.show()