Optimizing a function I'm optimizing the sum rate of two users in a communication system.
The problem can be formulated as follows:
$\begin{array}{l}
\mathop {\max }\limits_{\alpha ,\rho } {\rm{    }}\left( {1 + \frac{{\left( {1 - \rho } \right)\alpha P{g_1}}}{{\left( {1 - \rho  + \mu } \right){N_0}}}} \right)\left( {1 + \frac{{\left( {1 - \alpha } \right)P{g_2}}}{{\alpha P{g_2} + \left( {1 + \mu } \right){N_0}}} + \frac{{\rho \eta P{g_1}{g_3}}}{{\left( {1 + \mu } \right){N_0}}}} \right)\\
\text{subject to}\\
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \frac{{\left( {1 - \rho } \right)\left( {1 - \alpha } \right)P{g_1}}}{{\left( {1 - \rho } \right)\alpha P{g_1} + \left( {1 - \rho  + \mu } \right){N_0}}} \ge {T_2} \;\;\;\;\;\;\;\;\;\;\;\; (C1) \\
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \frac{{\left( {1 - \rho } \right)\alpha P{g_1}}}{{\left( {1 - \rho  + \mu } \right){N_0}}} \ge {T_1} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; (C2) \\
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\ \frac{{\left( {1 - \alpha } \right)P{g_2}}}{{\alpha P{g_2} + \left( {1 + \mu } \right){N_0}}} + \frac{{\rho \eta P{g_1}{g_3}}}{{\left( {1 + \mu } \right){N_0}}} \ge {T_2} \;\;\;\;\;\;\; (C3)\\
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\ 0 \le \alpha  \le 0.5 \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; (C4)\\
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\ 0 \le \rho  \le 1 \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; (C5)
\end{array}$
where $P$ is transmitted power, $N_0$ is noise power, ${g_i} = {\left| {{h_i}} \right|^2}$ with $h_i$ is the channel gain, $\mu$, $T_1$, $T_2$ are constants.
Without the constraints (C1)-(C3), we can prove that the objective function is maximized at $\alpha = 0$ or 0.5. Then, at $\alpha = 0$, $\rho$ will be 1. In addition, when $\alpha = 0.5$, the objective is a quadratic function, and thus, we can also find the solution analytically.
However, with the constraints (C1)-(C3), $\alpha$ and $\rho$ are coupled with each other in new constraints. This makes the problem much more difficult, and we are struggling to solve it.
Do you have any suggestions?
Thank you very much.
 A: One method is to use evolutionary algorithms such as Genetic Algorithm (GA). Which can be implemented in MATLAB or if the speed matters in C++.
Be careful about two points:


*

*Evolutionary algorithms are nondeterministic based on probability with no guarantee that the opimal point is a global optimum if the problem is non-convex (roughly blind search). However, the output is tried to be as best as possible.

*When using GA online, the initial population of the current step can be the final optimized population of the last step (warm start).
In the MATLAB implementation, you should code your variables into an array:
$$X=\begin{bmatrix}\alpha&\rho\end{bmatrix}^T$$
Here is a template example of setting a GA problem in MATLAB:
clear all
clc

tic
problem.nvars=2; % number of the variables
problem.Aineq=[]; % if there is linear inequality
problem.Bineq=[]; % if there is linear inequality
problem.Aeq=[]; % if there is linear equality
problem.Beq=[]; % if there is linear equality
problem.lb=[0,0]; % Variables lower bound
problem.ub=[0.5 1]; % Variables upper bound
problem.options = optimoptions(@ga);

problem.fitnessfcn= @my_cost; % the cost function
problem.nonlcon = @my_constraints % the constraint function
x_best = ga(problem);

toc

my_cost is a function written in file my_cost.m, it receives an $X=\begin{bmatrix}\alpha&\rho\end{bmatrix}^T$ and gives you the cost to be minimized. If you look for a maximization problem, simply negate the cost.
my_constraints is a function which checks the constraints and reject them if violate any rule. It is written in my_constraints.m files. For the implementation details see here.
The name of cost function and constraint function are arbitrary but the file name and the function name should match with each other.
There are also examples here.

I have written the C++ method. Just double check the formulation.
command to compile:
g++ -O3 -std=c++11 -pthread main.cpp

code:
#include <string>
#include <vector>
#include "genetic.hpp" // download from https://github.com/Arash-codedev/openGA
#include <fstream>
#include <sstream>
#include <iomanip>

double P=1.0; // ???
double T1=1.0; // ???
double T2=1.0; // ???
double g1=1.0; // ???
double g2=1.0; // ???
double g3=1.0; // ???
double N0=0.01; // ???
double mu=0.1; // ???
double etta=0.99; // ???

struct CommunicationOptions
{
    double alpha; // to be optimized
    double rho;   // to be optimized

    std::string to_string() const
    {
        std::ostringstream out;
        out<<"{";
        out<<"alpha: "<<alpha;
        out<<", ";
        out<<"rho: "<<rho;
        out<<"}";
        return out.str();
    }
};

struct MiddleCost
{
    // This is where the results of simulation
    // is stored but not yet finalized.
    double parenthesis1;
    double parenthesis2;
};

typedef EA::Genetic<CommunicationOptions,MiddleCost> GA_Type;
typedef EA::GenerationType<CommunicationOptions,MiddleCost> Generation_Type;

void init_genes(CommunicationOptions& p,const std::function<double(void)> &rand)
{
    // initialization:
    p.alpha=0.5*rand();
    p.rho=1.0*rand();
}

bool eval_genes(
    const CommunicationOptions& p,
    MiddleCost &c)
{
    // constexpr double pi=3.141592653589793238;
    // c.cost=10*double(p.x.size());
    // for(unsigned long i=0;i<p.x.size();i++)
    //  c.cost+=p.x[i]*p.x[i]-10.0*cos(2.0*pi*p.x[i]);
    double C1_left=(1.0-p.rho)*(1.0-p.alpha)*P*g1/((1.0-p.rho)*p.alpha*P*g1+(1.0-p.rho+mu)*N0);
    if(!(C1_left>=T2)) // violation of C1
        return false;
    double C2_left=(1.0-p.rho)*p.alpha*P*g1/((1.0-p.rho+mu)*N0);
    if(!(C2_left>=T1)) // violation of C2
        return false;
    double C3_left=(1.0-p.rho)*p.alpha*P*g2/(p.alpha*P*g2+(1.0+mu)*N0)+p.rho*etta*P*g1*g3/((1.0+mu)*N0);
    if(!(C3_left>=T2)) // violation of C3
        return false;
    c.parenthesis1=1.0+(1.0-p.rho)*p.alpha*P*g1/((1.0-p.rho+mu)*N0);
    c.parenthesis2=1.0+(1.0-p.rho)*p.alpha*P*g2/(p.alpha*P*g2+(1.0+mu)*N0)+p.rho*etta*P*g1*g3/((1.0+mu)*N0);
    return true;
}

CommunicationOptions mutate(
    const CommunicationOptions& X_base,
    const std::function<double(void)> &rand,
    double shrink_scale)
{
    // perform mutation:

    CommunicationOptions X_new;
    bool out_of_range;
    do{
        out_of_range=false;
        X_new=X_base;
        X_new.alpha+=(rand()*shrink_scale)*(rand()-rand());
        X_new.rho+=(rand()*shrink_scale)*(rand()-rand());
        if(X_new.alpha>0.5||X_new.alpha<0.0)
            out_of_range=true;
        if(X_new.rho>1.0||X_new.rho<0.0)
            out_of_range=true;
    } while(out_of_range);
    return X_new;
}

CommunicationOptions crossover(
    const CommunicationOptions& X1,
    const CommunicationOptions& X2,
    const std::function<double(void)> &rand)
{
    CommunicationOptions X_new;
    double r;
    r=rand();
    X_new.alpha=r*X1.alpha+(1.0-r)*X2.alpha;
    r=rand();
    X_new.rho=r*X1.rho+(1.0-r)*X2.rho;
    return X_new;
}

double calculate_SO_total_fitness(const GA_Type::thisChromosomeType &X)
{
    // finalize the cost
    double parenthesis1=X.middle_costs.parenthesis1;
    double parenthesis2=X.middle_costs.parenthesis2;
    double final_cost=parenthesis1*parenthesis2;
    return -final_cost;
}
std::ofstream output_file;

void SO_report_generation(
    int generation_number,
    const EA::GenerationType<CommunicationOptions,MiddleCost> &last_generation,
    const CommunicationOptions& best_genes)
{
    std::cout
        <<"Generation ["<<generation_number<<"], "
        <<"Best="<<last_generation.best_total_cost<<", "
        <<"Average="<<last_generation.average_cost<<", "
        <<"Best genes=("<<best_genes.to_string()<<")"<<", "
        <<"Exe_time="<<last_generation.exe_time
        <<std::endl;

    output_file
        <<generation_number<<"\t"
        <<last_generation.average_cost<<"\t"
        <<last_generation.best_total_cost<<"\t"
        <<best_genes.alpha<<"\t"
        <<best_genes.rho<<"\t"
        <<"\n";
}

int main()
{
    output_file.open("./output.txt");
    output_file
        <<"step"<<"\t"
        <<"cost_avg"<<"\t"
        <<"cost_best"<<"\t"
        <<"alpha"<<"\t"
        <<"rho"<<"\t"
        <<std::endl;

    EA::Chronometer timer;
    timer.tic();

    GA_Type ga_obj;
    ga_obj.problem_mode=EA::GA_MODE::SOGA;
    ga_obj.multi_threading=true;
    ga_obj.dynamic_threading=false;
    ga_obj.idle_delay_us=0; // switch between threads quickly
    ga_obj.verbose=false;
    ga_obj.population=10000; // 50-400 must be enough
    ga_obj.generation_max=1000;
    ga_obj.calculate_SO_total_fitness=calculate_SO_total_fitness;
    ga_obj.init_genes=init_genes;
    ga_obj.eval_genes=eval_genes;
    ga_obj.mutate=mutate;
    ga_obj.crossover=crossover;
    ga_obj.SO_report_generation=SO_report_generation;
    ga_obj.best_stall_max=40;
    ga_obj.average_stall_max=20;
    ga_obj.tol_stall_best=1e-6;
    ga_obj.tol_stall_average=1e-6;
    ga_obj.elite_count=10;
    ga_obj.crossover_fraction=0.7;
    ga_obj.mutation_rate=0.1;
    ga_obj.solve();

    std::cout<<"The problem is optimized in "<<timer.toc()<<" seconds."<<std::endl;

    output_file.close();
    return 0;
}

