该问题是一个标准的非线性非凸问题，因此该问题类别的任何求解器都适合解决该问题。

例如，以下代码在 MATLAB Toolbox YALMIP（免责声明，由我开发）中实现该问题，并使用局部非线性求解器 ipopt 解决该问题。

N = 3;
K = 3;
h = rand(N,K);
R = rand(K,1);
sigma = 1;
beta  = 1;
Pmax = 10;
Constraints = [];
P = sdpvar(N,K,'full');
C = sdpvar(N,K,'full');
sdpvar t
Constraints = [P>=0, C>=0, sum(sum(P)) <= Pmax,sum(C,1)<=1]
for k = 1:K
   Constraints = [Constraints, t<=(beta/(N*R(k)))*sum(C(:,k).*log2(1+P(:,k).*h(:,k)/sigma))];
end    
solvesdp(Constraints,-t,sdpsettings('solver','ipopt'))

ipopt 是一个局部求解器，所以你在这里只能期望一个局部最优解。但是，看起来您的问题实际上相当简单，因为解决方案通常似乎是全局最优的。您可以通过使用 YALMIP 中提供的全局求解器全局求解问题来看到这一点。它实现了一种基于线性松弛下界和非线性局部求解器上限的 b&b 策略。当然，它只适用于小问题（除非你有大量的空闲时间......）

solvesdp(Constraints,-t,sdpsettings('solver','bmibnb'))

最后，我想补充一点，您的问题结构实际上允许您消除 c 变量。可以看出，c的最优选择是使每一列中的所有元素都等于0，除了对数项向量中最大项对应的元素。因此，从概念上讲，您可以解决

Constraints = [P>=0, sum(sum(P))<= Pmax]
for k = 1:K
    Constraints = [Constraints, t<=(beta/(N*R(k)))*max(log2(1+P(:,k).*h(:,k)/sigma))];
end

当您在 YALMIP 中对其建模时，这会导致混合整数非线性非凸程序（引入二进制变量来模拟 max 算子的非凸使用），并且它似乎比直接的非线性方法效率低。