Hey guys,
I'm trying to implement a nonlinear constrained optimization with C++
using NLopt to compute one nash equilibrium of an n-person game. The
implementation is based on the function "npg" of the matlab optimization
toolbox, which solves an n-person finite non-co-operative game
(DOI-Link: http://dx.doi.org/10.1109/ICM2CS.2009.5397970). It uses SQP
based method to get a solution of the minimization problem.
The attached code represents my implementation, strictly oriented on the
matlab implementation referred above. Each call of the optimization
function returns "nlopt roundoff-limited" exception.
Does anybody know this exception or can figure out the fault? I'm not
very familiar with the mathematical background of different optimization
formulations. Maybe the calls of add_*_constraint or the formulation of
the objective function causes the exception.
You will need the Eigen library to compile the source code. It is
available at http://bitbucket.org/eigen/eigen/get/3.2.0.tar.bz2 .
Thanks in advance,
Tobias Schmidt
/*
* npg.cpp
*
* The function npg solves an n-person finite non-co-operative game to
* compute one sample Nash Equilibrium. It uses an optimization formulation
* of n-person non-co-operative games as described in the adjoining paper
* "An Optimization Formulation to Compute Nash Equilibrium in finite
* Games" presented by the author.
*
* The inputs to the function are:
* a) M : The row vector containing number of strategies of each player.
* b) U : The matrix containing the pay-offs to the players at various
* pure strategy combinations.
*
* The outputs are:
* a) A : The matrix whose columns are mixed strategies to players at Nash
* equilibrium.
* b) payoff : The row vector containing the pay-offs to the players at Nash
* equilibrium.
* c) iterations : Number of iterations performed by the optimization
* method.
* d) err : The absolute error in the objective function of the minimization
* problem.
*
* For theory of the method the adjoining paper may be consulted. Here an
* example is given for explanantion. Consider a 3-person game where each
* player has 2 strategies each. So M = [2 2 2]. Suppose the pay-offs at the
* pure strategy combinations (1-1-1), (1-1-2) and so on, as described by
* the ranking in the theory, are given as the matrix U =
* [2,7.5,0;3,.2,.3;6,3.6,1.5;2,3.5,5;0,3.2,9;2,3.2,5;0,2,3.2;2.1,0,1]. Then
* after supplying M and U call [A,payoff,iterations,err] = npg(M,U).
*
* The method is capable of giving one sample Nash equilibrium out of
* probably many present in a given game. The screenshot showing GUI has
* been developed on the code using it as dll and VB.Net. The GUI software
* may be made available on request.
*
* Any error in the code may be reported to bhaskerchatter...@gmail.com. Any
* suggestion/comment is greatly appreciated.
*
* Author: tschmidt
* See: n-person game for Optimization Toolbox on Matlab
*/
#include <string.h>
#include <cmath>
#include <iostream>
#include <nlopt.hpp>
#include "Eigen/Core"
using namespace Eigen;
typedef struct
{
const int s, p;
const MatrixXi& I;
const VectorXd& Us;
unsigned int myfun_calls;
} myfun_data_t;
typedef struct
{
const int s, p;
const VectorXi& pay;
const MatrixXi& I;
const MatrixXd& U;
} confun_data_t;
typedef struct
{
const MatrixXd& Aeq;
const VectorXd& beq;
const int s;
} linconfun_data_t;
double myfun(unsigned int n, const double* x, double* grad, void* f_data)
{
myfun_data_t* data = (myfun_data_t*) f_data;
double res = 0, prod = 1;
int i;
unsigned int j, k;
++data->myfun_calls;
for (k = 0; k < n - data->s; ++k)
res += x[data->s + k];
for (i = 0; i < data->p; ++i)
{
for (j = 0; j < n - data->s; ++j)
{
prod *= x[data->I(i,j)-1];
}
res -= data->Us(i) * prod;
prod = 1;
}
return res;
}
void confun(unsigned int m, double* result, unsigned int n, const double* x,
double* grad, void* f_data)
{
confun_data_t* data = (confun_data_t*) (f_data);
// Result vector
VectorXd C = VectorXd::Zero(m);
int i, j;
unsigned int t, k;
double add, prd;
for (i = 0; i < data->s; ++i)
{
C(i) = -x[data->pay(i)-1];
for (t = 0; t < n - data->s; ++t)
{
add = 0;
for (j = 0; j < data->p; ++j)
{
prd = 1;
for (k = 0; k < n - data->s; ++k)
if (i+1 == data->I(j,k))
prd *= data->U(j,k);
else
prd *= x[data->I(j,k)-1];
if (data->I(j,t) != i+1)
prd = 0;
add += prd;
}
C(i) += add;
}
}
// return ...
memcpy(result, C.data(), m * sizeof(double));
}
void linconfun(unsigned int m, double* result, unsigned int n, const double* x,
double* grad, void* f_data)
{
linconfun_data_t* data = (linconfun_data_t*) (f_data);
unsigned int i,j;
for (i = 0; i < n - data->s; ++i)
{
result[i] = -data->beq(i);
for (j = 0; j < m; ++j)
result[i] += data->Aeq(i,j) * x[j];
}
}
int npg(RowVectorXi& iM, MatrixXd& iU, MatrixXd& oA, RowVectorXd& oPayoff,
int* oIterations, double* oErr) throw(const std::string)
{
#if 1 // editor folding of npg matrices initialization code
int i = 0, j = 0, k = 0, p = 1;
double V = 1;
int n = iM.cols();
int s = iM.sum();
oA = MatrixXd::Zero(iM.maxCoeff(), n);
oPayoff = RowVectorXd::Zero(n);
for(i = 0; i < n; ++i)
p = p * iM(i);
if(p != iU.rows() || n != iU.cols())
throw std::string("Error: Dimension mismatch!"); //
ERR_DIM_MISMATCH;
RowVectorXd P = RowVectorXd::Zero(n);
RowVectorXd N = RowVectorXd::Zero(n);
P(n-1) = 1;
for(i = n-2; i >= 0; --i)
P(i) = P(i+1) * iM(i+1);
for(i = 0; i < n; ++i)
N(i) = p / P(i);
VectorXd x0 = VectorXd::Zero(s + n);
for(i = 0; i < n; ++i)
for(j=0; j < iM(i); ++j, ++k)
x0(k) = 1.0 / iM(i);
VectorXd Us = iU.rowwise().sum();
for(i = 0; i < n; ++i)
V = V * std::pow((1.0 / iM(i)), (double) iM(i));
x0.tail(n) << V * (iU.colwise().sum().transpose());
MatrixXd Aeq = MatrixXd::Zero(n, s+n);
int cnt = 0;
for (i = 0; i < n; ++i)
{
if (i != 0)
cnt += iM(i);
for (j = 0; j < s; ++j)
if ((j+1) <= iM.head(i+1).sum() && (j+1) > cnt)
Aeq(i,j) = 1;
}
VectorXd beq = VectorXd::Ones(n);
MatrixXi I = MatrixXi::Ones(p,n);
int counter = 0, count = 0;
for (i = 0; i < n; ++i)
for (j = 0; j < N(i); ++j)
{
++counter;
if (i != 0 && counter > iM.head(i+1).sum())
counter -= iM(i);
for (k = 0; k < P(i); ++k)
I(count++) = counter;
}
VectorXd lb = VectorXd::Zero(s+n);
VectorXd ub = VectorXd::Ones(s+n);
VectorXi pay = VectorXi::Zero(s);
counter = 0;
for (i = 0; i < n; ++i)
for (j = 0; j < iM(i); ++j)
pay(counter++) = i + 1 + s;
lb.tail(n).setConstant(-HUGE_VAL);
ub.tail(n).setConstant(HUGE_VAL);
#endif
// Construct optimization problem
nlopt::opt opt(nlopt::LD_SLSQP, s+n);
myfun_data_t fun_data = {s, p, I, Us, 0};
opt.set_min_objective(&myfun, &fun_data);
confun_data_t confun_data = {s, p, pay, I, iU};
opt.add_inequality_mconstraint(&confun, (void*) &confun_data,
std::vector<double>(s+n,1e-6));
linconfun_data_t linconfun_data = {Aeq, beq, s};
opt.add_equality_mconstraint(&linconfun, (void*) &linconfun_data,
std::vector<double>(s+n,1e-6));
std::vector<double> lbv(lb.data(), lb.data() + lb.rows());
std::vector<double> ubv(ub.data(), ub.data() + ub.rows());
opt.set_lower_bounds(lbv);
opt.set_upper_bounds(ubv);
opt.set_xtol_rel(1e-6);
opt.set_xtol_abs(1e-6);
opt.set_ftol_rel(1e-6);
opt.set_ftol_abs(1e-6);
opt.set_maxeval(1e4);
opt.set_maxtime(10.0);
// Call optimization function and map output
std::vector<double> x(x0.data(), x0.data() + x0.rows());
try
{
opt.optimize(x);
} catch (nlopt::roundoff_limited err)
{
throw std::string("Halted because roundoff errors limited
progress! ") + err.what();
} catch (nlopt::forced_stop err)
{
throw std::string("Halted because of a forced termination! ") +
err.what();
} catch (std::invalid_argument err )
{
throw std::string("Invalid arguments. ") + err.what();
} catch (std::bad_alloc err)
{
throw std::string("Ran out of memory. ") + err.what();
} catch (std::runtime_error err)
{
throw std::string("Generic failure. ") + err.what();
}
count = 0;
for (i = 0; i < n; ++i)
{
for (j = 0; j > iM(i); ++j)
oA(j,i) = std::abs(x[count++]);
oPayoff(i) = x[s+i];
}
*oErr = std::abs(opt.last_optimum_value());
*oIterations = fun_data.myfun_calls;
return EXIT_SUCCESS;
}
int main(int argc, char *argv[]) {
RowVectorXi m(2);
m << 3, 3;
MatrixXd u(9,2);
u << 0.1324, -0.1324,
0.1238, -0.1238,
0.1155, -0.1155,
0.1470, -0.1470,
0.1384, -0.1384,
0.1301, -0.1301,
0.1610, -0.1610,
0.1524, -0.1524,
0.1442, -0.1442;
MatrixXd a;
RowVectorXd payoff;
int iter = 0;
double err = 0.0;
try
{
npg(m, u, a, payoff, &iter, &err);
} catch(const std::string message)
{
std::cout << message << std::endl;
return EXIT_FAILURE;
}
std::cout << "Iterations: " << iter << "\nFunction value: " << err <<
std::endl;
return EXIT_SUCCESS;
}
_______________________________________________
NLopt-discuss mailing list
NLopt-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/nlopt-discuss
