Introduction
The update of the question invalidates my original approach, which depended on the problem being a simple special case. The current problem consists of a linear objective function with linear constraints, to be maximized over the integers. Theoretically, one can solve the problem with LinearProgramming. In the OP's setup, the variables betha form a (rank-2) matrix instead of a rank-1 vector needed to use LinearProgramming. One can flatten the problem out and use integer linear programming. Alternatively, one can do a brute-force search, which turns out to be faster in this case.
Brute force is a little faster
Define the cost in terms of a memoized coefficient function costCoeff, and it turns out an exhaustive search of the possible inputs is at least twice as fast as LinearProgramming. See further down for some discussion of the use of listability in the code for costCoeff etc. For instance, Times automatically multiplies corresponding parts of lists and matrices. Such properties allows us to rewrite the OP's code more efficiently in terms of the whole matrices and vectors instead of their parts.
constQ = 2; (* not given -- made up *)
P = 5; Ns = 7; G = 3;
SeedRandom[1];
usrK = RandomInteger[{2, 5}, G];
gh = RandomReal[{1, 3}, {G, Ns}];
costCoeff[P_?NumericQ, Ns_?NumericQ, gh_?(MatrixQ[#, NumericQ] &),
usrK_?(VectorQ[#, IntegerQ] &), G_?NumericQ] :=
costCoeff[P, Ns, gh, usrK, G] = -Exp[usrK/(P gh)] (usrK) / Log[2] ExpIntegralEi[-(usrK/(P gh))]
cost[P_?NumericQ, Ns_?NumericQ, gh_?(MatrixQ[#, NumericQ] &),
usrK_?(VectorQ[#, IntegerQ] &), G_?NumericQ, betha_] :=
Total[costCoeff[P, Ns, gh, usrK, G] betha, 2];
We use UnitStep[x - constQ] to implement the OP's qosConstr constraints. The values of UnitStep have be 1 on all rows for the constraints to be met; if one is not met, the minimum of the values will be 0. This criterion is used to Pick out of all the 0-1 matrices with exactly one 1 in each column, those that satisfy these constraints.
domain = Transpose[Tuples[Permutations[Append[ConstantArray[0, G - 1], 1]], 7], {1, 3, 2}];
domain2 =
Pick[
domain,
Min /@ Transpose @
UnitStep[Total[Transpose[domain, {3, 1, 2}] Log2[1 + (P/Ns) gh], {2}] - constQ],
1];
betamax = domain2[[
Last @ Ordering[cost[P, Ns, gh, usrK, G, Transpose[domain2, {3, 1, 2}]]]
]]
cost[P, Ns, gh, usrK, G, betamax]
(*
{{0, 0, 1, 0, 0, 1, 1}, {1, 0, 0, 0, 1, 0, 0}, {0, 1, 0, 1, 0, 0, 0}}
41.9077
*)
Linear programming is conceptually more natural
LinearProgramming minimizes its objective function. To get a maximum we can take the negative of the OP's cost function. Below, costCoeff represents the matrix of coefficients. To maximize the cost function, I used -Flatten @ costCoeff[P, Ns, gh, usrK, G] for the objective function. Each linear constraints is specified in two parts, a list of coefficients of the variables and a pair {b, s} representing a boundary value b and a code s -- -1, 0, or 1 -- for the inequality as described in the documentation for LinearProgramming. The second-to-last argument specifies the bounds on each variable to be 0 and 1.
solFlat = LinearProgramming[
-Flatten @ costCoeff[P, Ns, gh, usrK, G], (* objective function *)
Table[Flatten @ Transpose @ (* one 1 in each column *)
ReplacePart[ConstantArray[0., {Ns, G}], i -> ConstantArray[1, G]], {i, Ns}] ~Join~
MapIndexed[ (* fourth constraint coeffs *)
Flatten @ ReplacePart[ ConstantArray[0., {G, Ns}], #2 -> #1] &,
Log2[1 + (P/Ns) gh]],
ConstantArray[{1, -1}, Ns] ~Join~ (* constraint inequalities *)
ConstantArray[{constQ, 1}, G],
ConstantArray[{0, 1}, G*Ns], (* 0-1 constraint *)
Integers
];
sol = Partition[solFlat, Ns] (* reshape into G x Ns matrix *)
cost[P, Ns, gh, usrK, G, sol] (* maximum cost *)
(*
{{0, 0, 1, 0, 0, 1, 1}, {1, 0, 0, 0, 1, 0, 0}, {0, 1, 0, 1, 0, 0, 0}}
41.9077
*)
Remarks and explanation
Objective function - listability and efficiency
To use Mathematica efficiently with matrices, vectors, and other arrays, it pays to learn how theListable attribute (or more precisely, internal vectorization) and threading work with the numeric functions, including Plus and Times, as well as the importance of packed arrays (see What is a Mathematica packed array?). Here are a few ways to compute the OP's cost function, the OP's way costOP, the way from the original answer costME2, and the way above costFlat. They are equivalent.
costOP = Sum[-Exp[usrK[[g]]/(P gh[[g, n]])] (usrK[[g]] betha[g, n]) /
Log[2] ExpIntegralEi[-(usrK[[g]]/(P gh[[g, n]]))], {g, 1, G}, {n, 1, Ns}];
costME2 = Total[costCoeff[P, Ns, gh, usrK, G] Array[betha, {G, Ns}], 2];
costFlat = Flatten@costCoeff[P, Ns, gh, usrK, G] . Flatten@Array[betha, {G, Ns}];
costOP == costME2 == costFlat
(* True *)
Remark: In computing the objective function, instead of using Sum and Part, one should use the listability of the built-in numeric functions, which Mathematica handles much more efficiently. This would be an important point if the objective function were to be evaluated many times, as it would if NMinimize were used. In the above costCoeff is evaluated only once. If we turn them into functions, we can compare the timings:
ClearAll[costOP, costME2, costFlat];
costOP[betha_] := Sum[-Exp[usrK[[g]]/(P gh[[g, n]])] (usrK[[g]] betha[[g, n]])/
Log[2] ExpIntegralEi[-(usrK[[g]]/(P gh[[g, n]]))], {g, 1, G}, {n, 1, Ns}];
costFlat[betha_] := Flatten@costCoeff[P, Ns, gh, usrK, G] . Flatten@betha;
costME2[betha_] := Total[costCoeff[P, Ns, gh, usrK, G] betha, 2];
SeedRandom[1];
testOP = Table[costOP[RandomInteger[1, {G, Ns}]], {1000}]; // AbsoluteTiming
SeedRandom[1];
testFlat = Table[costFlat[RandomInteger[1, {G, Ns}]], {1000}]; // AbsoluteTiming
SeedRandom[1];
testME2 = costME2[Transpose[RandomInteger[1, {1000, G, Ns}], {3, 1, 2}]]; // AbsoluteTiming
(*
{0.372749, Null}
{0.130880, Null}
{0.000734, Null} <-- Wow. Compare with the brute-force method above.
*)
testOP == testFlat == testME2
(* True *)
Constraints
The constraints are implemented in three pieces.
First, the lists
Table[Flatten @ Transpose @ ReplacePart[ConstantArray[0., {Ns, G}], i -> {1, 1, 1}], {i, Ns}]
together with
ConstantArray[{1, -1}, Ns]
specify the column sums should be at most 1.
Second, the lists
MapIndexed[
Flatten @ ReplacePart[ConstantArray[0., {G, Ns}], #2 -> #1] &,
Log2[1 + (P/Ns) gh]]
together with
ConstantArray[{constQ, 1}, G]
specify the row sums with the coefficients from the corresponding row of the matrix Log2[1 + (P/Ns) gh] should be at most constQ.
Third, the argument
ConstantArray[{0, 1}, G*Ns], (* 0-1 constraint *)
specifies that each variable should be between 0 and 1 (inclusive), and the final argument Integers specifies that the variables are to be integers.
Timing
The AbsoluteTiming of the brute-force method is 0.001875 sec.
The AbsoluteTiming of the LinearProgramming code is 0.004235, which is more than twice as long as the brute-force method.
Original answer
First, your objective function is linear in your variables betha, which form a matrix in your setup.
Second, your constraints imply that an input is a 0-1 matrix with exactly one entry being 1 in each column. Edit note: the OP later added a constraint that makes the following not apply to the current question.
The fastest way to maximize the function is to pick among each column of coefficients, the position of the maximum coefficient and set the corresponding variable to 1 (and set the others in the column to zero, of course).
Transpose @
Map[ReplacePart[ConstantArray[0, G], # -> 1] &,
Last /@ Ordering /@ Transpose @ costCoeff[P, Ns, gh, usrK, G]
]
cost[P, Ns, gh, usrK, G, %]
(*
{{0, 1, 1, 0, 1, 1, 1}, {1, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}}
42.5
*)
A check
A simple check is to evaluate the function on all possible inputs and take the max.
domain = Transpose[
Tuples[Permutations[Append[ConstantArray[0, G - 1], 1]], 7], {1, 3, 2}];
domain[[
Last @ Ordering[cost[P, Ns, gh, usrK, G, Transpose[domain, {3, 1, 2}]]]
]]
(*
{{0, 1, 1, 0, 1, 1, 1}, {1, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}}
42.5
*)
