// (c) 2023-2024 Fair Isaac Corporation
using System;
using Optimizer;
using Optimizer.Objects;
using static Optimizer.Objects.Utils;
using static Optimizer.Objects.LinExpression;
namespace XpressExamples
{
/// <summary>QCQP problem (linear objective, convex quadratic constraints)</summary>
/// <remarks>
/// 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.
/// </remarks>
class Catenary
{
static readonly int N = 100; /// Number of chainlinks
static readonly int L = 1; /// Difference in x-coordinates of endlinks
static readonly double H = 2.0 * L / N; /// Length of each link
public static void Main(string[] args)
{
using (XpressProblem prob = new XpressProblem())
{
// Variables: x-coordinates
Variable[] x = prob.AddVariables(N + 1)
.WithName("x({0})")
.WithLB(XPRS.MINUSINFINITY)
.ToArray();
// Variables: y-coordinates
Variable[] y = prob.AddVariables(N + 1)
.WithName("y({0})")
.WithLB(XPRS.MINUSINFINITY)
.ToArray();
// Bounds: positions of endpoints
// Left anchor
x[0].Fix(0);
y[0].Fix(0);
// Right anchor
x[N].Fix(L);
y[N].Fix(0);
/// Objective: Minimise the potential energy
/// sum(j in 1..N) (y(j-1)+y(j))/2
LinExpression obj = LinExpression.Create();
for (int i = 1; i <= N; i++)
{
obj.AddTerm(0.5 * y[i - 1]).AddTerm(0.5 * y[i]);
}
prob.SetObjective(obj, ObjSense.Minimize);
/// Constraint: positions of chainlinks
/// forall(j in 1..N) (x(j)-x(j-1))^2+(y(j)-y(j-1))^2 <= H^2
for (int i = 1; i <= N; i++)
{
prob.AddConstraint(Sum(LinExpression.Create().AddTerm(x[i]).AddTerm(-x[i - 1]).Square(),
LinExpression.Create().AddTerm(y[i]).AddTerm(-y[i - 1]).Square())
.Leq(H * H).SetName($"Link_{i}"));
}
/// Solve the problem and write the solution
prob.LpOptimize("");
Console.WriteLine($"Solution: {prob.ObjVal}");
for (int i = 0; i <= N; i++)
{
Console.WriteLine($"{i}: {x[i].GetSolution()}, {y[i].GetSolution()}");
}
}
}
}
}
