OK, in my understanding the optimization should be executable with a gradient-based algorithm if I set *grad* to a constant value instead of a calculated one. A short try also fails with the nlopt roundoff- limited exception after the first optimization step.
I found out that the constant value, which I previously set as "gradient" was too small. A bigger one causes the program to run successfully! Hence, it led me implementing the gradient calculations based on the difference quotients.
Unfortunately, the optimization converges with an x vector which violates the constraints, but the "last_optimize_result()" says "FTOL_REACHED"!
Does either the way I define or add the functions cause the violations? Regards, Tobias
/* * 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 [email protected]. Any * suggestion/comment is greatly appreciated. * * Author: tschmidt * See: n-person game for Optimization Toolbox on Matlab */ #include <cstring> #include <math.h> #include <iostream> #include <cstdio> #include <nlopt.hpp> #include <fpstatus.h> #include <time.h> #include "Eigen/Core" using namespace Eigen; typedef struct { const int s, p; const MatrixXi& I; const VectorXd& Us; unsigned int myfun_calls; // double* pre_x; // double* pre_f; } myfun_data_t; typedef struct { const int p; const VectorXi& pay; const MatrixXi& I; const MatrixXd& U; // double* pre_x; // double* pre_f; } confun_m_data_t; typedef struct { const MatrixXi& Aeq; const VectorXi& beq; // double* pre_x; // double* pre_f; } linconfun_m_data_t; //typedef struct { // int x_off; // int x_end; // double c; //} linconfun_data_t; void printMatD(const MatrixXd &mat) { for (int i = 0; i < mat.rows(); ++i) { for (int j = 0; j < mat.cols(); ++j) std::printf("%f\t", mat(i, j)); std::printf("\n"); } } void printMatI(const MatrixXi &mat) { for (int i = 0; i < mat.rows(); ++i) { for (int j = 0; j < mat.cols(); ++j) std::printf("%d\t", mat(i, j)); std::printf("\n"); } } const double d = 1e-3; // delta to calculate difference quotient 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; unsigned int i, j; ++data->myfun_calls; if(grad){ std::cout << data->myfun_calls << ": myfun("; for (i = 0; i < n; ++i) std::printf("%.6f%s", x[i], i + 1 < n ? ", " : ") = "); } for (i = 0; i < n - data->s; ++i) res += x[data->s + i]; 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; } if(grad) std::cout << res << std::endl; if (grad) { double f1, *x1; x1 = (double*) std::malloc(n * sizeof(double)); std::memcpy(x1, x, n * sizeof(double)); for (i = 0; i < n; ++i) { x1[i] += d; f1 = myfun(n, x1, NULL, data); --data->myfun_calls; // correction of functions call counter grad[i] = (f1 - res) / d; x1[i] -= d; // static gradient // grad[j] = d; // gradient based on last vectors // grad[j] = (x[j] != data->pre_x[j]) ? (res - data->pre_f[j]) / (x[j] // - data->pre_x[j]) : 0; // data->pre_f[j] = res; // data->pre_x[j] = x[j]; } std::free(x1); } return res; } // Non-linear inequality constraint function void confun_m(unsigned int m, double* result, unsigned int n, const double* x, double* grad, void* f_data) { confun_m_data_t* data = (confun_m_data_t*) (f_data); unsigned int i, j, t, k; double add, prd; for (i = 0; i < m; ++i) { result[i] = -x[data->pay(i) - 1]; for (t = 0; t < n - m; ++t) { add = 0; for (j = 0; j < data->p; ++j) { prd = 1; for (k = 0; k < data->I.cols(); ++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; } result[i] += add; } // convert Constraint h(x) <= 0 to -h(x) >= 0 // result[i] = -result[i]; } // correct gradient if (grad) { double f1, *x1; x1 = (double*) std::malloc(n * sizeof(double)); std::memcpy(x1, x, n * sizeof(double)); for (i = 0; i < m; ++i) for (j = 0; j < n; ++j) { x1[j] += d; confun_m(1, &f1, n, x1, NULL, data); grad[i * n + j] = (f1 - result[i]) / d; x1[j] -= d; } std::free(x1); } // if (grad) { // for (k = 0; k < m; ++k) { // for (t = 0; t < n; ++t) { //// grad[k * n + t] = d; // grad[k * n + t] = (x[t] != data->pre_x[t]) ? (result[k] - data->pre_f[k * n // + t]) / (x[t] - data->pre_x[t]) : 0; // data->pre_f[k * n + t] = result[k]; // data->pre_x[t] = x[t]; // } // } // } } // Linear Equality constraint function void linconfun_m(unsigned int m, double* result, unsigned int n, const double* x, double* grad, void* f_data) { linconfun_m_data_t* data = (linconfun_m_data_t*) (f_data); unsigned int i, j; for (i = 0; i < m; ++i) { result[i] = -data->beq(i); for (j = 0; j < n; ++j) result[i] += data->Aeq(i, j) * x[j]; } // correct gradient if (grad) { double f1, *x1; x1 = (double*) std::malloc(n * sizeof(double)); std::memcpy(x1, x, n * sizeof(double)); for (i = 0; i < m; ++i) for (j = 0; j < n; ++j) { x1[j] += d; linconfun_m(1, &f1, n, x1, NULL, data); grad[i * n + j] = (f1 - result[i]) / d; x1[j] -= d; } std::free(x1); } // grad based on last values // if (grad) { // for (i = 0; i < m; ++i) { // for (j = 0; j < n; ++j) { //// grad[i * n + j] = d; // grad[i * n + j] = (x[j] != data->pre_x[j]) ? (result[i] - data->pre_f[i * n // + j]) / (x[j] - data->pre_x[j]) : 0; // data->pre_f[i * n + j] = result[i]; // data->pre_x[j] = x[j]; // } // } // } } //double linconfun(unsigned n, const double* x, double* grad, void* f_data) //{ // linconfun_data_t* data = (linconfun_data_t*) f_data; // double res = data->c; // unsigned int i; // for (i = data->x_off; i < data->x_end; ++i) // res -= x[i]; // // if (grad) { // double f1, *x1; // x1 = (double*) std::malloc(n * sizeof(double)); // std::memcpy(x1, x, n * sizeof(double)); // for (i = 0; i < n; ++i) { // x1[i] += d; // f1 = linconfun(n, x1, NULL, data); // grad[i] = (f1 - res) / d; // x1[i] -= d; // } // std::free(x1); // } // return res; //} nlopt::result 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()); // x0.tail(n) << (iU.colwise().sum().transpose()) / iU.rows(); // Use this line for COBYLA MatrixXi Aeq = MatrixXi::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; } VectorXi beq = VectorXi::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 #if 1 // nlopt optimization problem // Construct optimization problem nlopt::opt opt(nlopt::LD_SLSQP, s + n); // nlopt::opt opt(nlopt::LN_COBYLA, s+n); // double* fun_pre_x = (double*) std::calloc(s + n, sizeof(double)); // double* fun_pre_f = (double*) std::calloc(s + n, sizeof(double)); myfun_data_t fun_data = { s, p, I, Us, 0 /*, fun_pre_x, fun_pre_f */}; opt.set_min_objective(&myfun, &fun_data); // double* confun_pre_x = (double*) std::calloc(s + n, sizeof(double)); // double* confun_pre_f = (double*) std::calloc(s * (s + n), sizeof(double)); confun_m_data_t confun_m_data = { p, pay, I, iU /*, confun_pre_x, confun_pre_f */}; opt.add_inequality_mconstraint(&confun_m, (void*) &confun_m_data, std::vector<double>(s, 1e-10)); // double* linconfun_pre_x = (double*) std::calloc(s + n, sizeof(double)); // double* linconfun_pre_f = // (double*) std::calloc(n * (s + n), sizeof(double)); linconfun_m_data_t linconfun_m_data = { Aeq, beq /*, linconfun_pre_x, linconfun_pre_f */}; opt.add_equality_mconstraint(&linconfun_m, (void*) &linconfun_m_data, std::vector<double>(n, 1e-10)); 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-2); // opt.set_ftol_rel(1e-6); opt.set_ftol_abs(1e-6); // opt.set_maxeval(200); // opt.set_maxtime(.09); // Call optimization function and map output std::vector<double> x(x0.rows()); try { x = opt.optimize(std::vector<double>(x0.data(), x0.data() + x0.rows())); } 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(); } #endif count = 0; for (i = 0; i < n; ++i) { for (j = 0; j < iM(i); ++j) oA(j, i) = std::fabs(x[count++]); oPayoff(i) = x[s + i]; } *oErr = std::fabs(opt.last_optimum_value()); *oIterations = fun_data.myfun_calls; return opt.last_optimize_result(); } 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; // RowVectorXi m(3); // m << 3, 3, 3; // MatrixXd u(27, 3); // u << -28.0266, 56.5718, -28.5452, -28.0266, 56.6301, -28.6035, -28.0266, 56.6839, -28.6573, -28.0622, 56.5726, -28.5103, -28.0622, 56.6307, -28.5685, -28.0622, 56.6845, -28.6222, -28.0952, 56.5734, -28.4782, -28.0952, 56.6316, -28.5363, -28.0952, 56.6852, -28.5900, -27.9673, 56.5125, -28.5452, -27.9673, 56.5708, -28.6035, -27.9673, 56.6245, -28.6573, -28.0029, 56.5132, -28.5103, -28.0029, 56.5714, -28.5685, -28.0029, 56.6251, -28.6222, -28.0358, 56.5140, -28.4782, -28.0358, 56.5721, -28.5363, -28.0358, 56.6258, -28.5900, -27.9126, 56.4578, -28.5452, -27.9126, 56.5160, -28.6035, -27.9126, 56.5698, -28.6573, -27.9481, 56.4584, -28.5103, -27.9481, 56.5166, -28.5685, -27.9481, 56.5703, -28.6222, -27.9809, 56.4592, -28.4782, -27.9809, 56.5173, -28.5363, -27.9809, 56.5709, -28.5900; // u << -27.5591, 56.3438, -28.7846, -27.5591, 56.4057, -28.8465, -27.5591, 56.4634, -28.9042, -27.5974, 56.3449, -28.7476, -27.5974, 56.4067, -28.8094, -27.5974, 56.4644, -28.8670, -27.6331, 56.3462, -28.7131, -27.6331, 56.4080, -28.7748, -27.6331, 56.4656, -28.8324, -27.5041, 56.2887, -28.7846, -27.5041, 56.3506, -28.8465, -27.5041, 56.4083, -28.9042, -27.5422, 56.2898, -28.7476, -27.5422, 56.3516, -28.8094, -27.5422, 56.4093, -28.8670, -27.5780, 56.2910, -28.7131, -27.5780, 56.3528, -28.7748, -27.5780, 56.4104, -28.8324, -27.4540, 56.2386, -28.7846, -27.4540, 56.3005, -28.8465, -27.4540, 56.3582, -28.9042, -27.4921, 56.2396, -28.7476, -27.4921, 56.3014, -28.8094, -27.4921, 56.3591, -28.8670, -27.5277, 56.2408, -28.7131, -27.5277, 56.3025, -28.7748, -27.5277, 56.3601, -28.8324; MatrixXd a; RowVectorXd payoff; int iter = 0; double err = 0.0; clock_t begin, end; nlopt::result opt_return; try { std::cout << "Running npg..." << std::endl; if ((begin = clock()) == -1) throw std::string("Error getting clock time"); opt_return = npg(m, u, a, payoff, &iter, &err); if ((end = clock()) == -1) throw std::string("Error getting clock time"); } catch (const std::string message) { std::cerr << message << std::endl; return EXIT_FAILURE; } std::cout << "Optimization returned with "; switch (opt_return) { case nlopt::SUCCESS: std::cout << "SUCCESS" << std::endl; break; case nlopt::STOPVAL_REACHED: std::cout << "STOPVAL_REACHED" << std::endl; break; case nlopt::FTOL_REACHED: std::cout << "FTOL_REACHED" << std::endl; break; case nlopt::XTOL_REACHED: std::cout << "XTOL_REACHED" << std::endl; break; case nlopt::MAXEVAL_REACHED: std::cout << "MAXEVAL_REACHED" << std::endl; break; case nlopt::MAXTIME_REACHED: std::cout << "MAXTIME_REACHED" << std::endl; break; default: std::cout << "an error: " << opt_return << std::endl; } std::cout << "Objective function calls: " << iter << "\nAbsolute error: " << err << "\nTime: " << (end - begin) / (double) CLOCKS_PER_SEC << std::endl; std::cout << "A:" << std::endl; printMatD(a); std::cout << "Payoff:" << std::endl; printMatD(payoff); return EXIT_SUCCESS; }
_______________________________________________ NLopt-discuss mailing list [email protected] http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/nlopt-discuss
