FICO Xpress Optimization Examples Repository
 FICO Optimization Community FICO Xpress Optimization Home

Catenary - Solving a QCQP

Description
This model finds the shape of a hanging chain by minimizing its potential energy.

 xbcatenac.zip [download all files]

Source Files
By clicking on a file name, a preview is opened at the bottom of this page.
 xbcatena.c [download]

xbcatena.c

/********************************************************
BCL Example Problems
====================

file xbcatena.c

QCQP problem (linear objective, convex quadratic constraints)
Based on AMPL model catenary.mod
(Source: http://www.orfe.princeton.edu/~rvdb/ampl/nlmodels/ )

This model finds the shape of a hanging chain by
minimizing its potential energy.

(c) 2008-2023 Fair Isaac Corporation
author: S.Heipcke, June 2008, rev. Mar. 2011
********************************************************/

#include <stdio.h>
#include "xprb.h"
#include "xprs.h"

#define N 100	               /* Number of chainlinks */
#define L 1	               /* Difference in x-coordinates of endlinks */

double H = 2.0*L/N;	       /* Length of each link */

int main(int argc, char **argv)
{
int i;
XPRBprob prob;
XPRBvar x[N+1], y[N+1];       /* x-/y-coordinates of endpoints of chainlinks */
XPRBctr cobj, c;

prob=XPRBnewprob("catenary");           /* Initialize a new problem in BCL */

/**** VARIABLES ****/
for(i=0;i<=N;i++)
x[i] = XPRBnewvar(prob, XPRB_PL, XPRBnewname("x(%d)",i), -XPRB_INFINITY,
XPRB_INFINITY);
for(i=0;i<=N;i++)
y[i] = XPRBnewvar(prob, XPRB_PL, XPRBnewname("y(%d)",i), -XPRB_INFINITY,
XPRB_INFINITY);
/* Bounds: positions of endpoints */
/* Left anchor */
XPRBfixvar(x[0],0);  XPRBfixvar(y[0],0);
/* Right anchor */
XPRBfixvar(x[N],L);  XPRBfixvar(y[N],0);

/****OBJECTIVE****/
/* Minimise the potential energy:  sum(j in 1..N) (y(j-1)+y(j))/2  */
cobj = XPRBnewctr(prob, "Obj", XPRB_N);
for(i=1;i<=N;i++)
{
XPRBaddterm(cobj, y[i-1], 0.5);
XPRBaddterm(cobj, y[i], 0.5);
}
XPRBsetobj(prob, cobj);                 /* Set the objective function */

/**** CONSTRAINTS ****/
/* Positions of chainlinks:
forall(j in 1..N) (x(j)-x(j-1))^2+(y(j)-y(j-1))^2 <= H^2  */
for(i=1;i<=N;i++)
{
c = XPRBnewctr(prob, XPRBnewname("Link_%d",i), XPRB_L);
XPRBaddqterm(c, x[i], x[i], 1);
XPRBaddqterm(c, x[i], x[i-1], -2);
XPRBaddqterm(c, x[i-1], x[i-1], 1);
XPRBaddqterm(c, y[i], y[i], 1);
XPRBaddqterm(c, y[i], y[i-1], -2);
XPRBaddqterm(c, y[i-1], y[i-1], 1);
XPRBaddterm(c, NULL, H*H);
}

/****SOLVING + OUTPUT****/
XPRBsetsense(prob, XPRB_MINIM);         /* Choose the sense of optimization */

/* Problem printing and matrix output: */
/*
XPRBprintprob(prob);
XPRBexportprob(prob, XPRB_MPS, "catenary");
XPRBexportprob(prob, XPRB_LP, "catenary");
*/

/* Start values:
for(i=0;i<=N;i++) XPRBsetinitval(x[j], j*L/N);
for(i=0;i<=N;i++) XPRBsetinitval(y[j], 0);
*/
/* Disable convexity check */
XPRSsetintcontrol(XPRBgetXPRSprob(prob), XPRS_IFCHECKCONVEXITY, 0);

XPRBlpoptimize(prob,"");               /* Solve the problem */

printf("Solution: %g\n", XPRBgetobjval(prob));
for(i=0;i<=N;i++)
printf(" %d: %g, %g\n", i, XPRBgetsol(x[i]), XPRBgetsol(y[i]));

return 0;
}