# Quadratic programming for special equation issues

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;

{
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>");
}

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:

using System;

{
class Program
{
const double alpha = 0.2;
const int k = 10;

static void Main(string[] args)
{
const double tSumMax = 0.04;
const int loops = 10;
const int t1Loops = 1000;
const int t2Loops = 1000;
double t1min = -0.002;
double t1max = 0.02;
double t2min = -0.002;
double t2max = 0.025;
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++)
{
double t1 = t1min + i1 * dt1;

for (int i2 = 0; i2 < t2Loops; i2++)
{
double t2 = t2min + i2 * dt2;

if (t1 + t2 <= tSumMax)
{
double val = E(t1, t2);

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>");
}

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(rho, t1, t2);
}

return r / k;
}

static double R(double rho, double t1, double t2)
{
double r = t1 * (10 - rho) + t2 * (10 - rho) * rho/10
+ (t1 + t2)*(10 - rho) * Math.Pow(rho/10, 20);

return r;
}

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(rho, t1, t2) - _M;

v += s * s;
}

return v / k;
}
}
}


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\}$.

• Great code. Thank you so much – John Nov 5 '14 at 3:51
• Thank Axel Kemper. What is function that you shown? For the old function, you are correct. However, I add some new term in the old function. The result is not change. Could you see my new function and let me know if it is correct? – John Nov 5 '14 at 14:57
• Yes, I used the "old" function. The new function is quite different from the old one and much more involved. Sorry, but I won't have the time to dig through your code. – Axel Kemper Nov 5 '14 at 15:41

Why would you want to solve it using sequential quadratic programming? The problem is a convex quadratic program. Hence, any standard convex QP solver will solve it trivially (alternatively, you can solve it easily by hand, just derive the quadratic expression, check stationary point. If not feasible, look at the boundary of the simple polytope you are optimizing over)

The following piece of code sets up and solves the problem in MATLAB using the toolbox YALMIP (disclaimer, developed by me).

sdpvar tau1 tau2
k = 10;
M = 0;
for rho = 0:k-1
R = tau1*(10-rho)+tau2*(10-rho)*(rho/10) - (tau1+tau2)*(10-rho)*(rho/10)^20;
M = M + R;
end
M = M/k;
V = 0;
for rho = 0:k-1
R = tau1*(10-rho)+tau2*(10-rho)*(rho/10) - (tau1+tau2)*(10-rho)*(rho/10)^20;
V = V + (R-M)^2;
end
V = V/k;
alpha = .2;
E = M - alpha*V % Linear - Convex quadratic

Constraints = [-.002 <= tau1 <= 0.02, -.002 <= tau2 <= .025,tau1+tau2 <= .04];
% Maximize E, i.e., minimize -E
optimize(Constraints,-E)
% ans =

% 0.0200
% 0.0200

• Sorry. I cannot run your code by matlab. It must setup something, right? I wwant to use QP to resolve my problem to compare with the original paper – John Nov 5 '14 at 9:14
• You mean you have MATLAB but cannot run the code? To run it, you must install YALMIP (google it) – Johan Löfberg Nov 5 '14 at 9:29
• I run it. Please see my update function. Is it possible to find the my new update formula. I also update matlab code for find $\Omega$ at dropbox.com/s/qs2egatyaiqjtog/robustsoliton.m?dl=0 To run it please run by [omega, ddeg, cumddist] = robustsoliton(500,0.02,0.01) – John Nov 5 '14 at 10:00
• It is missing a function/variable idealsoliton. However, the question is still the same, since $R$ is linear in $\tau$. You haven't changed anything really, from an optimization point of view. – Johan Löfberg Nov 5 '14 at 10:08
• That is exactly what I am saying. Applying find on a symbolic variable ($\tau$) leads to an extremely complicated set of constraints to represent, so it is not allowed. However, if I simply set d_max=k, it computed value([tau1 tau2]) as [-0.0062 0.0262] – Johan Löfberg Nov 5 '14 at 10:47