My problem is how to find $\tau_1$ and $\tau_2$ s.t maximize the objective function is $$E=M-\alpha V$$ subject to $$-0.0062\le\tau_1\le0.499$$ $$-0.479\le\tau_2\le0.0262$$ $$\tau_1+\tau_2\le0.02$$ where$$M=\frac{1}{k}\sum^{k-1}_{\rho=0}R(\rho)$$ $$V=\frac{1}{k}\sum^{k-1}_{\rho=0}(R(\rho)-M)^2$$ $\alpha$ is learning rate and sets equals $0.27$
Let see my example in which $k=500$: $$R(\rho)=(1+0.1) \times [\tau_1\times f(\rho,1)+\tau_2\times f(\rho,2)-(\tau_1+\tau_2)\times f(\rho,221)+\frac {k-\rho}{1+0.1}\ln{\frac{k-\rho}{k}}+\sum _{d=1}^{221} \Omega_d\times f(\rho,d)]$$ In which $f(\rho,d)=(k-\rho)\times d \times (\frac{\rho}{k})^{d-1}$
$\Omega_d$ is a element in probability distribution $\Omega$. I done to make it
My question is that can we use Sequential quadratic programming to find the solution of objective function? Please give me some suggestion. I would like to thank Mr. Axel Kemper for your suggestion. Based on the your suggestion, I modify my code as bellow. However, my result is not similar the optimal result of reference paper $(0.006,0.026)$. My result is $(0.002,-0.002)$. Have some problem in here?
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace QuadraticProgramming
{
class Program
{
const double alpha = 0.27;
const int k = 500;
public static double[] omega;
public static int d_max = k;
static void Main(string[] args)
{
double c = 0.02;
double delta = 0.01;
omega = RSD(k, c, delta);
d_max = 221;
double tSumMax = omega[d_max];
const int loops = 100;
const int t1Loops = 100;
const int t2Loops = 100;
double t1min = -0.0062;
double t1max = 0.499;
double t2min = -0.479;
double t2max = 0.0262;
double t1opt = (t1max + t1min) / 2;
double t2opt = (t2max + t2min) / 2;
double maxVal = E(t1opt, t2opt);
for (int i = 0; i < loops; i++)
{
double dt1 = (t1max - t1min) / t1Loops;
double dt2 = (t2max - t2min) / t2Loops;
// scan search area for maximum
for (int i1 = 0; i1 < t1Loops; i1++)
{
// Console.WriteLine("{0}", i1);
double t1 = t1min + i1 * dt1;
for (int i2 = 0; i2 < t2Loops; i2++)
{
double t2 = t2min + i2 * dt2;
if (t1 + t2 <= tSumMax && (t1min <= t1) && (t1 <= t1max) && (t2min <= t2) && (t2 <= t2max))
{
//Console.WriteLine("dmax={0}, omega={1} ", d_max,omega[d_max]);
double val = E(t1, t2);
//Console.WriteLine("dmax={0},maxVal={1} ", d_max, maxVal);
if (val > maxVal)
{
maxVal = val;
t1opt = t1;
t2opt = t2;
Console.WriteLine("{0} {1} {2}", t1, t2, val);
}
}
}
}
// shrink search area
t1min = Math.Max(t1min, t1opt - dt1);
t1max = Math.Min(t1max, t1opt + dt1);
t2min = Math.Max(t2min, t2opt - dt2);
t2max = Math.Min(t2max, t2opt + dt2);
}
Console.WriteLine("<press any key to quit>");
Console.ReadKey();
}
static double E(double t1, double t2)
{
return M(t1, t2) - alpha * V(t1, t2);
}
static double M(double t1, double t2)
{
double r = 0;
for (int rho = 0; rho < k; rho++)
{
r += R(k, d_max, omega, rho, t1, t2);
}
return r / k;
}
static double R(int k, int d_max, double[] omega, double rho, double t1, double t2)
{
double sumterm = 0;
for (int d = 1; d <= d_max; d++)
{
sumterm += omega[d - 1] * f_d(k, rho, d);
}
double r = t1 * f_d(k, rho, 1) + t2 * f_d(k, rho, 2)
- (t1 + t2) * f_d(k, rho, d_max) + sumterm + (k - rho) / (double)(1 + 0.1) * Math.Log((k - rho) / (double)k);
return (1 + 0.1) * r;
}
static double f_d(int k, double rho, double d)
{
double f = (k - rho) * d * Math.Pow(rho / (double)k, (d - 1));//d=0=>d=1
return f;
}
static double V(double t1, double t2)
{
double _M = M(t1, t2);
double v = 0;
for (int rho = 0; rho < k; rho++)
{
double s = R(k, d_max, omega, rho, t1, t2) - _M;
v += s * s;
}
return v / k;
}
// Create Omega distribution of 1 to k with parameters c and delt
static double[] RSD(int k, double c, double delta)
{
double S = c * Math.Log(k / (double)delta) * Math.Sqrt((double)k);
double[] rsdMatrix = new double[k];
double[] rho = new double[k];
double[] tau = new double[k];
double sumRT = 0;
int kS = 0;
int i = 0;
for (i = 0; i < k; i++)
{
tau[i] = 0;
}
rho[0] = 1 / (double)k;
for (i = 1; i < k; i++)
{
rho[i] = 1 / (double)(i * (i + 1));//i=0<=>d=1=>d-1=i
}
kS = (int)Math.Round(k / (double)S);
if (kS > k)
kS = k;
for (i = 0; i < kS - 1; i++)
{
tau[i] = S / (double)k * 1 / (double)(i + 1);
}
tau[kS] = S / (double)k * Math.Log(S / (double)delta);
for (i = kS + 1; i < k; i++)
{
tau[i] = 0;
}
for (i = 0; i < k; i++)
{
sumRT += rho[i] + tau[i];
}
for (i = 0; i < k; i++)
{
rsdMatrix[i] = (rho[i] + tau[i]) / sumRT;
}
return rsdMatrix;
}
}
}
I have implemented a brute-force search procedure in
C#to find optimal values for $\tau_1$ and $\tau_2$.Within the allowed domain, all value pairs for $\tau_1$ and $\tau_2$ are tried along a grid and the best pair remembered.
The search domain is shrunk after every loop to improve accuracy. Function $R(\rho,\tau_1,\tau_2)$ should be rewritten to cope with values $k \ne 10$.
After some seconds of runtime, the procedure comes up with $0.02$ for both variables.
My
C#code:To actually see how well-behaved the function is, I created a contour plot:
The plot looks very regular indeed. It actually confirms that the maximum is at $\{0.02, 0.02\}$.