Catenary: Determine chain shape

# Example for a QCQP problem (linear objective, convex quadratic constraints).
#
# This model finds the shape of a hanging chain by
# minimizing its potential energy.
#
# Based on AMPL model catenary.mod:
#    (Source: https://vanderbei.princeton.edu/ampl/nlmodels/catenary/)
#
# (C) 2018-2025 Fair Isaac Corporation

import xpress as xp

N = 100                       # Number of chainlinks.
L = 1                         # Difference in x-coordinates of endlinks.
H = 2*L/N                     # Length of each link.

RN = range(N+1)

p = xp.problem()

x = p.addVariables(N+1, lb=-xp.infinity, name="x")
y = p.addVariables(N+1, lb=-xp.infinity, name="y")

# Objective: minimise the potential energy.
p.setObjective(xp.Sum((y[j-1] + y[j]) / 2 for j in range(1, N+1)))

# Bounds: positions of endpoints.
# Left anchor.
p.addConstraint(x[0] == 0)
p.addConstraint(y[0] == 0)
# Right anchor.
p.addConstraint(x[N] == L)
p.addConstraint(y[N] == 0)

# Constraints: positions of chainlinks.
p.addConstraint((x[j] - x[j-1])**2 + (y[j] - y[j-1])**2 <= H**2
                for j in range(1, N+1))

# Uncomment to export the matrix file.
# p.writeProb('catenary.mat', 'l')

p.optimize()

print("Solution: ", p.attributes.objval)
xsol, ysol = p.getSolution(x), p.getSolution(y)
for j in RN:
    print("{0:10.5f} {1:10.5f}".format(xsol[j],
                                       ysol[j]))