Open JianqiangDing opened 1 year ago
yuanchenyang
commented
1 year ago Can you check that the inputs from picos into mosek remain the same between runs? If your problem has multiple solutions, any difference in the inputs to mosek may result in a different solution. Another way to force an unique solution is to add an objective to the problem.
JianqiangDing
commented
1 year ago Hi, yuanchen,
I just finish the tests you suggested, I add an objective as sum of all coefficient squares name gamma, and write out the problem into files before and after calling prob.solve, you may check these these files, and there are some small numerical differences between them, like
I also tried to write out the problem by mosek, but it failed (error code 4000) because of
MSK_RES_ERR_INVALID_FILE_FORMAT_FOR_SYM_MAT (4000)
The file format does not support a problem with symmetric matrix variables.Here is the code, at the end of this code are two results, and you may notices coefficients of these two polynomials are different in some degree, especially the constant term.
import picos
import sympy as sp
from SumOfSquares import SOSProblem, Basis, prod
from sympy import Matrix, sin
from picos import norm
import numpy as np
def polynomial(name, variables, deg, hom=False):
"""Returns a (possibly homogeneous) degree DEG polynomial in VARIABLES,
with a variable (a sympy symbol) named using NAME for each coefficient. Used
in Sum-of-Squares relaxations for polynomial optimization.
"""
variables = sorted(variables, key=str) # use lex order
basis = Basis.from_degree(len(variables), deg, hom=hom)
coeffs = sp.symbols(f"{name}_:{len(basis)}")
p = sum(
coeff * prod(var ** power for var, power in zip(variables, monom))
for monom, coeff in zip(basis, coeffs)
)
return p, coeffs, basis
if __name__ == "__main__":
x, y, u = sp.symbols("x y u")
f = Matrix(
[
[x - 0.01 * (0.5 * x + 0.5 * y - 0.5 * x * y)],
[y + 0.01 * (-0.5 * y + 1 + sin(u))],
]
)
h = x ** 2 + y ** 2 - 1
h0 = x ** 2 + y ** 2 - 1.1
g = 10 * x ** 2 + 10 * (y - 0.5) ** 2 - 1
point = [0, -0.5]
epsilon = 1e-8
lambda_ = 1.01
degree = 6
ds = 2
prob = SOSProblem()
var_xy = [x, y]
V, VC, VB = polynomial("coe", var_xy, degree)
VC = [prob.sym_to_var(c) for c in VC]
gamma = 0
for c in VC:
gamma += c ** 2
v1 = V.subs([(x, f[0, 0]), (y, f[1, 0])])
VE = sp.integrate(v1, (u, -sp.pi, sp.pi)) / (2 * sp.pi)
V0 = VE.subs([(x, point[0]), (y, point[1])])
s0, coe0, _ = polynomial("coe0", var_xy, ds)
s1, coe1, _ = polynomial("coe1", var_xy, ds)
s2, coe2, _ = polynomial("coe2", var_xy, ds)
s3, coe3, _ = polynomial("coe3", var_xy, ds)
s4, coe4, _ = polynomial("coe4", var_xy, ds)
prob.add_sos_constraint(VE - lambda_ * V + s0 * h - s1 * g, var_xy)
prob.add_sos_constraint(-V + s2 * h0 - s3 * h, var_xy)
prob.add_sos_constraint(V0 - epsilon, var_xy)
prob.add_sos_constraint(s0, var_xy)
prob.add_sos_constraint(s1, var_xy)
prob.add_sos_constraint(s2, var_xy)
prob.add_sos_constraint(s3, var_xy)
prob.add_sos_constraint(s4, var_xy)
prob.set_objective(direction='min', expression=gamma)
prob.write_to_file('test_pre_result01.dat-s', writer='picos')
prob.solve(solver="mosek")
prob.write_to_file('test_after_result01.dat-s', writer='picos')
variables = sorted(var_xy, key=str) # use lex order
coef = [c.value for c in VC]
p = sum(
coeff * prod(var ** power for var, power in zip(variables, monom))
for monom, coeff in zip(VB, coef)
)
print(p)
import matplotlib.pyplot as plt
delta = 0.005
x = np.arange(-1.5, 1.5, delta)
y = np.arange(-1.5, 1.5, delta)
X, Y = np.meshgrid(x, y)
p = sp.simplify(p)
f = sp.lambdify(var_xy, p, "numpy")
Z = f(X, Y)
px = 1 / plt.rcParams["figure.dpi"]
width, height = 800, 800
fig, ax = plt.subplots(figsize=(width * px, height * px), layout="constrained")
CS = ax.contour(X, Y, Z, [0], colors="black")
f = sp.lambdify(var_xy, h, "numpy")
Z = f(X, Y)
CS = ax.contour(X, Y, Z, [0], colors="red")
f = sp.lambdify(var_xy, g, "numpy")
Z = f(X, Y)
CS = ax.contour(X, Y, Z, [0], colors="green")
ax.autoscale_view()
ax.axis("equal")
ax.axis("off")
plt.show()
# -9.78540636846771e-6 * x ** 6 - 6.54459698023813e-6 * x ** 5 * y - 5.49648030440343e-6 * x ** 5
# - 3.31036623986582e-6 * x ** 4 * y ** 2 + 4.99430964535994e-6 * x ** 4 * y - 6.78898583577452e-6 * x ** 4
# - 2.59906609471962e-6 * x ** 3 * y ** 3 - 3.42960814456312e-6 * x ** 3 * y ** 2
# - 2.41911575511092e-6 * x ** 3 * y - 1.39066673181138e-6 * x ** 3 - 8.44765011180873e-6 * x ** 2 * y ** 4
# + 3.81630450700703e-6 * x ** 2 * y ** 3 - 6.33074859235012e-6 * x ** 2 * y ** 2
# + 1.85576039236168e-6 * x ** 2 * y - 2.20676546027284e-5 * x ** 2 - 5.74773409048156e-6 * x * y ** 5
# - 5.94900095136526e-6 * x * y ** 4 - 2.16945782913507e-6 * x * y ** 3 - 6.83476880990713e-6 * x * y ** 2
# + 4.82373411338204e-9 * x * y + 1.77024516308475e-6 * x - 5.95443197187209e-6 * y ** 6
# + 3.94073215280344e-6 * y ** 5 - 3.82728277687708e-7 * y ** 4 - 1.57335448618497e-6 * y ** 3
# + 5.22005453462501e-7 * y ** 2 - 2.83448824424391e-7 * y - 1.89368929498165e-7
# -9.80079936955942e-6 * x ** 6 - 6.51463408446183e-6 * x ** 5 * y - 5.46705161393525e-6 * x ** 5
# - 3.32552055316688e-6 * x ** 4 * y ** 2 + 4.96683338395406e-6 * x ** 4 * y - 6.81459176944389e-6 * x ** 4
# - 2.60042011878491e-6 * x ** 3 * y ** 3 - 3.38991519093161e-6 * x ** 3 * y ** 2
# - 2.34912922720159e-6 * x ** 3 * y - 1.34615560524517e-6 * x ** 3 - 8.38689792279441e-6 * x ** 2 * y ** 4
# + 3.78287202437611e-6 * x ** 2 * y ** 3 - 6.38212953812153e-6 * x ** 2 * y ** 2
# + 1.89315055407056e-6 * x ** 2 * y - 2.2003050952949e-5 * x ** 2 - 5.76312748873915e-6 * x * y ** 5
# - 5.89388305555828e-6 * x * y ** 4 - 2.17212703069132e-6 * x * y ** 3 - 6.79204268662201e-6 * x * y ** 2
# - 1.4294930220735e-8 * x * y + 1.74598287887084e-6 * x - 6.07836513027998e-6 * y ** 6
# + 4.07426303039097e-6 * y ** 5 - 4.50033444220666e-7 * y ** 4 - 1.65486376120298e-6 * y ** 3
# + 6.60483794235288e-7 * y ** 2 - 2.99434193124444e-7 * y - 2.27658259624662e-7I looked into this again, seems like the non-determinism might be caused by some floating point issues in sympy, as a lot of your values are very close to zero. They don't happen for any of my test cases.
JianqiangDing
commented
1 year ago Hi, chenyang,
yea, I agree with you, your cases work perfectly and this may be caused by sympy.
I also checked YALMIP later, and the outputs look all same (according to the final plots), this may heavily rely on MATLAB's powerful symbolic computation. And I also tried setting picos numerical tolerance to be compatible with the backend solver as described in Numeric Tolerances, but it failed.
I suggest you may leave this issue open for now for anyone who may be facing the same problems, and we may fix this issue in the near future:-).
Hi, yuanchen,
I am using SumOfSquares to reproduce the results of a paper, the problem is when I finish the program locally, the solution may be varying each time when I run it, Is this a bug of SumOfSquares or just some numerical issue? In my own opinion, It is unreasonable this shall be a numerical issue because of you may get some solution deviate to the exact answer, but you shall not get different wrong answer each time without introducing randomness.
I only notice that you control the numerical tolerance manually at https://github.com/yuanchenyang/SumOfSquares.py/blob/5edf9ca7ebfc3478802d900ede4408bbad07b467/src/SumOfSquares/util.py#L35
and I changed this locally, I guess this probably not the core of the problem.
And if you change the variable from
ds=2tods=6, the results seems more stable each time, but according to the printed result, you can also known they are not all same each running.Any suggestion would helpful, thank you.
Here is the code