I have an optimization problem to solve using Quantlib. However, the documentation (or the lack of it) is confusing me no end.. Can someone please provide some pointers or code snippets which will at least help me in the right direction ?I have an array a[ a1,a2,a3,a4] . Lets say a1 = 0.1, a2 = 0.2, a3 = 0.3, a4 = 0.4. such that a1 + a2 + a3 + a4 = 1.0The values can be within +-10% of the initial values. So, bounds of a1 = 0.09-0.11, bounds of a2 = 0.18-0.22 and so on..I have another weights array w[ w1,w2,w3,w4].I just have to maximise the dot product (a1w1 + a2w2 + a3w3 + a4w4) subject to the condition that (a1 + a2 + a3 + a4 ) = 1.0.Any pointers please ?

You probably want to use something like GLPK (http://ftp.gnu.org/gnu/glpk/) rather than QuantLib for this type of linear programming problem.

I cannot use GLPK primarily because I need a less restrictive license. Can someone help me please Quantlib experts ?

Something like the code attached might be a start. It compiles and run on my machine, but I haven't checked at all that the results are correct.Sorry I didn't have the time to document it; if you can't make sense of it, ask away and I'll try to clarify.Also, most parameters are hard-coded for lack of time, but I hope it won't be hard for you to generalize.Luigi

Thank you for your response! I actually coded this up yesterday using Simplex optimizer in a manner very similar to what you have done. It indeed works.

Hello Luigi, and other experts,I coded up this solution; and ran it with ConjugateGradient, BFGS and Simplex optimizers. Unfortunately the results seem to be erratic, and are not consistent. They seem to be actually incorrect in case of CD/BFGS for certain inputs; and only true for Simplex for certain values of Lambda. It would be very beneficial if someone could let me know if I can work these out consistently with CG/BFGS, OR, a way to locate appropriate Lambda so that it works for all inputs with Simplex.Again to reiterate the problem : ======================I have an array a[ a1,a2,a3,a4] . Lets say a1 = 0.1, a2 = 0.2, a3 = 0.3, a4 = 0.4. such that a1 + a2 + a3 + a4 = 1.0The values can be within +-20% of the initial values. So, bounds of a1 = 0.09-0.11, bounds of a2 = 0.18-0.22 and so on..I have another weights array w[ w1,w2,w3,w4].I just have to minimise the dot product (a1w1 + a2w2 + a3w3 + a4w4) subject to the condition that (a1 + a2 + a3 + a4 ) = 1.0 (the initial input sum)My Program:=========#include <cmath>#include <vector>#include <iostream>#include <ql/math/optimization/constraint.hpp>#include <ql/math/optimization/costfunction.hpp>#include <ql/math/optimization/simplex.hpp>#include <ql/math/optimization/conjugategradient.hpp>#include <ql/math/optimization/bfgs.hpp> #include <ql/math/optimization/levenbergmarquardt.hpp>using namespace QuantLib;using namespace std;void maxAlpha (std::vector<double> &vec, std::vector<double> &alpha, double Pct);class AlphaCostFunction : public CostFunction { public : AlphaCostFunction(const Array &Arr, double iSum ) { alpha = Arr; InputSum = iSum; } double value ( const Array & x) const { double sum = 0; double res = DotProduct(x, alpha); for (unsigned int i = 0 ; i < alpha.size() ; i++) sum += x; res += (sum - InputSum ) * (sum - InputSum) * 100000.0; //I AM DOING THIS SO THAT THE RESULT SUMS TO INPUTSUM. return res; } Disposable <Array > values ( const Array & x) const { std::cout << "NEVER REALLY COMES HERE!" << std::endl; Array res (13); for (int m = 0 ; m < 13 ; m++) res[m]= x[m]; return res; } private: Array alpha; double InputSum;};class AlphaConstraint : public Constraint { public: class Impl : public Constraint::Impl { public: Impl(const Array &V, const double&Pct) { initVal = V; initPct = Pct; std::cout << "Lower Bounds: " << std::endl; for (unsigned int i = 0 ; i < initVal.size() ; i++) std::cout << (initVal - ((Pct/100) * initVal)) << " "; std::cout << endl; std::cout << "Upper Bounds: " << std::endl; for (unsigned int i = 0 ; i < initVal.size() ; i++) std::cout << (initVal + ((Pct/100) * initVal)) << " "; std::cout << endl; } bool test(const Array& x) const { bool ret = true; for (unsigned int i = 0 ; (i < x.size()) && ret ; i++) { ret = ret && ((initVal - (initPct * std::abs(initVal) / 100)) <= x) && (x <= (initVal + (initPct * std::abs(initVal) / 100))); // I AM DOING THIS TO IMPOSE BOUNDS } return ret; } private: Array initVal; double initPct; }; AlphaConstraint(const Array& initV, const double& initPct) : Constraint(boost::shared_ptr<Constraint::Impl>(new AlphaConstraint::Impl(initV, initPct))) {} };void maxAlpha (std::vector<double> &vec, std::vector<double> &alpha, double Pct) { double iSum = 0; Array xInp(vec.size()); Size maxIterations =20000; Size minStatIterations =100; double Epsilon =1e-12; EndCriteria myEndCrit ( maxIterations , minStatIterations, Epsilon, Epsilon, Epsilon); for (unsigned int i = 0 ; i < vec.size() ; i++) { xInp = vec; iSum += vec; } std::cout << "Input Array is " << xInp << std::endl; double pct = Pct; Array alphaArray(alpha.size()); for (unsigned int i = 0 ; i < alpha.size() ; i++) alphaArray = alpha; std::cout << "Alpha Array is " << alphaArray << std::endl; AlphaCostFunction myFunc(alphaArray, iSum) ; AlphaConstraint constraint(xInp, pct); Problem myProb ( myFunc , constraint , xInp); Array scaleArray(vec.size()); ConjugateGradient solver1; BFGS solver2; Simplex solver3(0.5); EndCriteria :: Type solvedCrit = solver1.minimize ( myProb , myEndCrit ); if (solvedCrit != EndCriteria::None) { std::cout << "ConjugateGradient found a solution! " << std::endl; std :: cout << " Criteria :" << solvedCrit << std :: endl ; std :: cout << " Root :" << myProb . currentValue () << std :: endl ; std :: cout << " Min F Value :" << myProb . functionValue () << std :: endl ; } solvedCrit = solver2.minimize ( myProb , myEndCrit ); if (solvedCrit != EndCriteria::None) { std::cout << "BFGS found a solution! " << std::endl; std :: cout << " Criteria :" << solvedCrit << std :: endl ; std :: cout << " Root :" << myProb . currentValue () << std :: endl ; std :: cout << " Min F Value :" << myProb . functionValue () << std :: endl ; } solvedCrit = solver3.minimize ( myProb , myEndCrit ); if (solvedCrit != EndCriteria::None) { std::cout << "Simplex found a solution! " << std::endl; std :: cout << " Criteria :" << solvedCrit << std :: endl ; std :: cout << " Root :" << myProb . currentValue () << std :: endl ; std :: cout << " Min F Value :" << myProb . functionValue () << std :: endl ; }}int main(){ std::vector<double> vec; vec.push_back(1.95); vec.push_back(2.28); vec.push_back(2.03); vec.push_back(2.50); vec.push_back(2.51); vec.push_back(2.32); vec.push_back(2.39); vec.push_back(2.88); vec.push_back(2.82); vec.push_back(3.59); vec.push_back(4.05); vec.push_back(4.69); vec.push_back(9.79); std::vector<double> alpha; alpha.push_back(1); alpha.push_back(2); alpha.push_back(3); alpha.push_back(4); alpha.push_back(5); alpha.push_back(6); alpha.push_back(7); alpha.push_back(8); alpha.push_back(9); alpha.push_back(10); alpha.push_back(11); alpha.push_back(12); alpha.push_back(13); double Pct = 20.0; maxAlpha(vec, alpha, Pct);}OUTPUT:=======Input Array is [ 1.95; 2.28; 2.03; 2.5; 2.51; 2.32; 2.39; 2.88; 2.82; 3.59; 4.05; 4.69; 9.79 ] Alpha Array is [ 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13 ] Lower Bounds: 1.56 1.824 1.624 2 2.008 1.856 1.912 2.304 2.256 2.872 3.24 3.752 7.832 Upper Bounds: 2.34 2.736 2.436 3 3.012 2.784 2.868 3.456 3.384 4.308 4.86 5.628 11.748 ConjugateGradient found a solution! Criteria :StationaryFunctionValue Root :[ 2.34; 2.60499; 2.28998; 2.69498; 2.63997; 2.38496; 2.38995; 2.81494; 2.68993; 3.39493; 3.78992; 4.36491; 9.3999 ] Min F Value :366.425 BFGS found a solution! Criteria :StationaryFunctionValue Root :[ 2.34; 2.60499; 2.28998; 2.69498; 2.63997; 2.38496; 2.38995; 2.81494; 2.68993; 3.39493; 3.78992; 4.36491; 9.3999 ] Min F Value :366.425 Simplex found a solution! Criteria :StationaryPoint Root :[ 2.34; 2.72501; 2.40623; 2.88884; 2.83798; 2.56188; 2.48377; 3.03099; 2.79727; 3.44675; 3.59011; 3.75594; 8.9351 ] Min F Value :358.123As the Alpha Array is in an ascending order, the first element of the result MUST clip to the first upper bound, which it does correctly. However, the last value must also clip to lower bound [7.832] which it doesn't. None of the solutions are optimal. Queries:-Can the program be further optimized so that all three give correct and optimized results ?-If not, for some values of Lambda, Simplex does give correct results. How to locate the correct Lambda for a group of input numbers which will always give correct results ?Any response will be lot appreciated.Thanks!

Hi SabanyI think you should use a linear programming optimiser to solve your problem. Good, free LP optimisers exist, e.g. GNU Linear Programming Tool Kit or lpsolve. A linear optimiser will be much faster than all gradient/simplex based optimisers for your kind of problem and it will produce reliable results.IMO your problem is that gradiend based optimiser can be used in conjunction with box constraints (e.g. a1 <= 0.5), but if you have a constraint like a1 + a2 + a3 + a4 = 1.0which defines a hyperplan, normal implementations of gradiend/simplex based optimiser often fail (nothing special with QL).

I'm not an expert, but might you try constraining the weights to add to 1 instead of adding a penalty for not being 1? (See the cost function in my first example, where you let the first three weights vary and get the last as the difference from 1.) This would also reduce the dimension of the problem, which might help.