//+------------------------------------------------------------------+
//|                 slsqp_optimization_with_analytical_gradients.mq5 |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#property version   "1.00"
#include<slsqp_article\slsqp.mqh>
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
//---custom constraint parameters
   vector aa = {2,-1};
   vector bb = {0, 1};
//---gradient configuration instance
   GradDiffOptions obj_options;
//---set upper and lower limits for optimization problem
   obj_options.bounds = matrix::Zeros(2,2);
   obj_options.bounds.Fill(double("inf"));
   obj_options.bounds[0,0]*=-1.;
   obj_options.bounds[1,0] = 1.e-6;
//---configure gradient calculation method
   obj_options.method = GRAD_POINT_CALLABLE;
//---objective function object instance
   CObjectiveFunc obj_func;
//---inequality constraints function object instance
   CIneqConstraints ineq_func;
//---initial guess
   vector initp = {1.234, 5.678};
//---pass extra arguments to objective function's constraints  
   ineq_func.SetInternalData(aa,bb);
//---configure gradient options for objective function and constraints
   obj_func.setParams(obj_options);
   ineq_func.setParams(2,2,obj_options);
//---initialize and verify configuration of objective function and corresponding constraints
   if(!obj_func.initialize(initp) || !ineq_func.initialize(initp))
      return;
//---optimizer instance 
   CSlsqp slsqp_minim;
//---configure the optimizer
   slsqp_minim.SetXtolRel(1.e-4);
   vector ineq_tols(2);
   ineq_tols.Fill(1.e-8);
   slsqp_minim.SetInequalityConstraints(ineq_func,ineq_tols);
//---run the optimizer
   OptimizeResult ro = slsqp_minim.Minimize(obj_func);
//---check optimizer return code
   Print("Minimization return code : ", ro.return_code);
//---print results of optimization
   if(ro.return_code>0)
     {
      Print("Solution ", ro.solution);
      Print("Minimum value of function ", ro.objective_result);
      Print("Gradient at the minimum ", ro.objective_gradient);
      Print("Number of iterations ", ro.niter);
      Print("Number of function evaluations ", ro.nfeval);
     }
  }
//+------------------------------------------------------------------+
//+------------------------------------------------------------------+
//|objective function object                                         |
//+------------------------------------------------------------------+
class CObjectiveFunc:public CFunctor
  {
private:
    vector obj_grad;
public:
                     CObjectiveFunc(void)
     {
      obj_grad = vector::Zeros(2);
     }
                    ~CObjectiveFunc(void)
     {
     }
   virtual double    orig_fun(vector& x) override
     {
      return sqrt(x[1]);
     }
   virtual vector    grad_fun(vector& x) override
     {
      obj_grad[0] = 0.0;
      obj_grad[1] = 0.5/sqrt(x[1]);
      return obj_grad;
     }
  };
//+------------------------------------------------------------------+
//| inequality constraints                                           |
//+------------------------------------------------------------------+
class CIneqConstraints: public CConstraints
  {
private:
   vector            a, b;
   vector            f_result;
   matrix            g_result;
public:
                     CIneqConstraints(void)
     {
      f_result = vector::Zeros(2);
      g_result = matrix::Zeros(2,2);
     }
                    ~CIneqConstraints(void)
     {
     }
   void              SetInternalData(vector& a_, vector& b_)
     {
      a = a_;
      b = b_;
     }
   virtual vector    orig_fun(vector& x) override
     {
      f_result[0] = (a[0]*x[0] + b[0]) * (a[0]*x[0] + b[0]) * (a[0]*x[0] + b[0]) - x[1];
      f_result[1] = (a[1]*x[0] + b[1]) * (a[1]*x[0] + b[1]) * (a[1]*x[0] + b[1]) - x[1];
      return f_result;
     }
   virtual matrix    orig_grad(vector& x) override
     {
      g_result[0,0] = 3 * a[0] * (a[0]*x[0] + b[0]) * (a[0]*x[0] + b[0]);
      g_result[0,1] = -1.;
      g_result[1,0] = 3 * a[1] * (a[1]*x[0] + b[1]) * (a[1]*x[0] + b[1]);
      g_result[1,1] = -1.;
      return g_result;
     }
  };