I'm trying to fit a polynomial of the third degree through a number of points. This could be a very simple problem when not constraining the derivative. I found some promising solutions using CVXPY to combine least squares with constraints but so far I'm not even able to solve the least-squares problem with CVXPY. I also posted this first on Computational Science but given the fact that it is more about the code Stack Overflow might be more appropriate. The code is below.
import cvxpy as cp
import numpy as np
# Problem data
x = np.array([100, 200, 300, 400, 500])
y = 2160 - 1.678571429 * x + 0.006785714 * x * x + 5.52692E-20 * x * x * x
print(x)
print(y)
# Constructing the problem
a = cp.Variable(1)
b = cp.Variable(1)
c = cp.Variable(1)
d = cp.Variable(1)
objective = cp.Minimize(cp.sum_squares(a + b * x + c * x * x + d * x * x * x - y))
prob = cp.Problem(objective)
# The optimal objective value is returned by `prob.solve()`.
result = prob.solve(solver=cp.SCS, eps=1e-5)
# The optimal value for x is stored in `x.value`.
print("a: ", a.value)
print("b: ", b.value)
print("c: ", c.value)
print("d: ", d.value)
print(a.value + b.value * x + c.value * x * x + d.value * x * x * x)
The result of all this is a = 831.67009518; b = 1.17905623; c = 0.00155167 and d = 2.2071452e-06. This gives the mediocre result of y_est = [967.29960654 1147.20547453 1384.63057033 1692.81776513 2085.00993011]. Given that it should be able to yield a perfect solution, the provided solution is not really satisfying to me.
Once this works, I should be able to add the constraints as follows.
# add constraint for derivative = 0 in x = 500
constraints = [b + 2 * c * 500 + 3 * d * 500 * 500 == 0]
prob = cp.Problem(objective, constraints)
I'm not bound to CVXPY, so alternative solutions are also welcome.
It seems you just have a problem with association and how numpy and cvxpy differ in what * means. For example, c * x * x is not the same as x * x * c. The former is of course (c * x) * x and the second * is a dot product and thus the expression is a scalar. The latter (x * x) * c is what you want, as it first does an element-wise multiply. After doing that you get a much more sensible solution :)
a: [2159.91670117]
b: [-1.67850696]
c: [0.00678545]
d: [-1.13389519e-12]
[2059.92053699 2095.63343178 2267.05537873 2574.18637106 3017.02640194]