This is not going to be an answer per se as it is not going to offer a solution to the question asked.
It will simply demonstrate reasons for thinking the question itself more thoroughly.
The function in the question is a special case of $f(\vec x)=\sum_{i=1}^{k}a_ix_i+b\sum_{i=1}^{k}x_ilog(\frac{x_i}{\sum_{j=1}^{k}x_j})$, where $a_i<0$ and $b>0$.
The first derivatives of such functions are given by $\frac{\partial f}{\partial x_i}=a_i+blog(\frac{x_i}{\sum_{j=1}^{k}{x_j}})$.
The Hessian matrix ($H$) of such functions has on its main diagonal $b(\frac{1}{x_i}-\frac{1}{\sum_{j=1}^{k}{x_j}})$ and off-diagonal elements are equal to $-\frac{b}{\sum_{j=1}^{k}{x_j}}$. Note, that $Det(H)\equiv0$.
Additionally, in the question-where $k=6$-there are $4$ equality constraints, $6$ inequality constraints and $6$ non-negativity constraints.
If we assume that the proportions of each class of constraints remain constant, as the size of the problem increases, then the following code produces random instances of problems similar to the one in the question:
BlockRandom[
Block[{vars, pattn, sum, logs, alphas, betas, obj, A1, b1, A2, lambdas, mus, ypsilons, eqcs, ineqs, nncs, l},
(* 'Nn' is the number of variables and 'Ss' is the symbol to use for the variables *)
With[{Nn = 3, Ss = "n"},
(* returns a list of variables *)
vars[n_, str_] := Table[Unique[str, Temporary], n];
(* returns a list of variables' patterns *)
pattn[vars_] := Pattern[#, Blank[]] & /@ vars;
With[{vs = vars[Nn, Ss]},
With[{vspattns = pattn[vs]},
(* returns the sum of its arguments *)
sum[vspattns] := Evaluate[Total[vs]];
(* returns a list of 'Log[argument/sum[arguments]]' *)
logs[vspattns] := Evaluate[Log[vs/sum[vs]]];
With[{sm = sum[vs], lgs = logs[vs]},
With[{mu = -140223.9568567524, s = 110397.30093092154, cst = 8920.922},
(* original constraints *)
alphas = {-78416.91961691371, -293404.0587042522, -36099.19105260982, -53730.85220261606, -261977.41085974447,-117715.30870437803};
betas = {8920.922, 8920.922, 8920.922, 8920.922, 8920.922, 8920.922};
(* override these definitions to define a different (non-random) objective function *)
alphas = RandomVariate[NormalDistribution[mu, s], Nn];
betas = ConstantArray[cst, Nn];
(* define the objective function *)
obj[Sequence @@ vspattns] := Evaluate[Total[alphas vs + betas vs lgs]]
];
(* 'neqs' equals the number of (random) equality constraints to use *)
With[{neqcs = Round[Nn 2./3]},
(* original constraints *)
A1 = {{1, 1, 1, 0, 0, 0}, {1, 2, 0, 0, 1, 0}, {0, 0, 4, 2, 2, 0}, {0, 0, 0, 0, 0, 2}};
b1 = {3, 4, 8, 2};
A2 = {{0, 1, 1, 1, 1, 1}, {1, 0, 1, 1, 1, 1}, {1, 1, 0, 1, 1, 1}, {1, 1, 1, 0, 1, 1}, {1, 1, 1, 1, 0, 1}, {1, 1, 1, 1, 1, 0}};
(* override these definitions to define different (non-random) linear constraints *)
A1 = RandomVariate[DiscreteUniformDistribution[{0, 4}], {neqcs, Nn}];
b1 = RandomVariate[DiscreteUniformDistribution[{3, 8}], neqcs];
A2 = ConstantArray[1, {Nn, Nn}] - IdentityMatrix[Nn];
(* define the multipliers for the equality constraints *)
lambdas = vars[neqcs, "λ"];
(* define the multipliers for the inequality constraints *)
mus = vars[Nn, "μ"];
(* define the multipliers for the non-negativity constraints *)
ypsilons = vars[Nn, "υ"];
Block[{lambdaspattns, muspattns, ypsilonspattns},
(* define patterns associated with the multipliers of the equality constraints *)
lambdaspattns = pattn[lambdas];
(* define patterns associated with the multipliers of the inequality constraints *)
muspattns = pattn[mus];
(* define patterns associated with the multipliers of the non-negativity constraints *)
ypsilonspattns = pattn[ypsilons];
Block[{eqcsvs, ineqsvs, nncsvs, lvs, eqcspattns, ineqspattns, nncspattns, lpattns},
(* define relevant variables for the equality constraints equations *)
eqcsvs = Flatten[{vs, lambdas}];
(* define relevant variables for the inequality constraints equations *)
ineqsvs = Flatten[{vs, mus}];
(* define relevant variables for the non-negativity constraints equations *)
nncsvs = Flatten[{vs, ypsilons}];
(* define relevant variables for the Lagrangian *)
lvs = Flatten[{vs, lambdas, mus, ypsilons}];
(* define patterns for the definition of the equality constraints equations *)
eqcspattns = Flatten[{vspattns, lambdaspattns}];
(* define patterns for the definition of the inequality constraints equations *)
ineqspattns = Flatten[{vspattns, muspattns}];
(* define patterns for the definition of the non-negativity constraints equations *)
nncspattns = Flatten[{vspattns, ypsilonspattns}];
(* define patterns for the definition of the Lagrangian *)
lpattns = Flatten[{vspattns, lambdaspattns, muspattns, ypsilonspattns}];
(* The following definitions assist in constructing the Lagrangian *)
(* definition of the equality constraints equations *)
eqcs[Sequence @@ eqcspattns] := Evaluate[Flatten[(A1.Transpose[{vs}] - Transpose[{b1}]) lambdas]];
(* definition of the inequality constraints equations *)
ineqs[Sequence @@ ineqspattns] := Evaluate[Flatten[(A2.Transpose[{vs}]) mus]];
(* definition of the non-negativity constraints equations *)
nncs[Sequence @@ nncspattns] := Evaluate[ypsilons vs];
(* definition of the Lagrangian *)
l[Sequence @@ lpattns] := Evaluate[obj[Sequence @@ vs] - Total[eqcs[Sequence @@ eqcsvs]] - Total[ineqs[Sequence @@ ineqsvs]] - Total[nncs[Sequence @@ nncsvs]]];
(* obtain relations conducive to a minimum *)
Block[{foceqs, eqeqs, csceqs1, csceqs2, sys},
(* first order conditions *)
foceqs = Thread[D[l @@ lvs, {vs}] == 0.] // Simplify;
(* equality constraints *)
eqeqs = Thread[0. == (eqcs @@ eqcsvs)/lambdas] // Simplify;
(* complementary slackness conditions wrt to inequality constraints *)
csceqs1 = Thread[0. == (ineqs @@ ineqsvs)] // Simplify;
(* complementary slackness conditions wrt to non-
negativity constraints *)
csceqs2 = Thread[0. == (nncs @@ nncsvs)];
(* all equations bundled together *)
sys = Flatten[{foceqs, eqeqs, csceqs1, csceqs2}];
sys
]
]
]
]
]
]
]
]
], RandomSeeding -> 123]
In the code, matrices A1, b1, A2 are used to reproduce the constraints; changing their relative position, can recreate the original problem (also please remember to change the number of variables Nn to 6 in that case.)
After evaluation, the code returns a list of equations that should hold true at a minimum (it does not return the inequalities that should hold true at the minimum namely, the inequality and non-negativity constraints of the original problem and the non-negativity constraints of the respective multipliers).
However, that system of equations-either for the original problem or of the random problems I tested-cannot be solved. On my machine, I get
"Solve::ratnz: Solve was unable to solve the system with inexact
coefficients. The answer was obtained by solving a corresponding exact
system and numericizing the result."
and then an empty list.
I do not think that this is an issue with Mathematica, rather it has to do with the problem itself; the vanishing Hessian is probably an indication that something might go wrong. Perhaps, using the Lagrangian is not the ideal way to manage this optimization problem.
ps. the relevant inequalities can be assembled using
{
Thread[(ineqs @@ ineqsvs)/mus >= 0],
Thread[(nncs @@ nncsvs)/ypsilons >= 0],
Thread[mus >= 0],
Thread[ypsilons >= 0]
}
before or after sys, in the code.
FindRootinstead ofSolve. $\endgroup$Minimizemight be able to use Lagrange multipliers natively. Am I right in understanding thatg1==g2==g3==g4==0? $\endgroup$Reduceon them returnsFalse. $\endgroup$