/********************************************************
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 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;
}
|
