The following program uses the Davidson-Fletcher-Powell (DFP) method to find the minimum of a function. This method is a quasi-Newton method. That is, the DFP algorithm is based on Newton's method but performs different calculations to update the guess refinements.
Click here to download a ZIP file containing the project files for this program.
The program prompts you to either use the predefined default input values or to enter the following:
1. The initial guesses for the optimum point for each variable..
2. The tolerance for each variable,
3. The tolerance for the function
4. The maximum number of iteration cycles
In case you choose the default input values, the program displays these values and proceeds to find the optimum point. In the case you select being prompted, the program displays the name of each input variable along with its default value. You can then either enter a new value or simply press Enter to use the default value. This approach allows you to quickly and efficiently change only a few input values if you so desire.
The program displays the following final results:
1. The coordinates of the minimum point.
2. The minimum function value.
3. The number of iterations
The current code finds the minimum for the following function:
f(x1,x2) = x1 - x2 + 2 * x1 ^ 2 + 2 * x1 * x2 + x2 ^ 2
Using, for each variable, an initial value of 0, initial step size of 0.1, minimum step size of 1e-7, and using a function tolerance of 1e-7. Here is the sample console screen:
Here is the listing for the main module. The module contains several test functions:
using System.Diagnostics;
using System.Data;
using System.Collections;
using Microsoft.VisualBasic;
using System.Collections.Generic;
using System;
namespace Optim_DFP1
{
sealed class Module1
{
static public void Main()
{
int nNumVars = 2;
double[] fX = new double[] { 0.5, 0.5 };
double[] fParam = new double[] { 0, 0 };
double[] fToler = new double[] { 0.00001, 0.00001 };
int nIter = 0;
int nMaxIter = 100;
double fEpsFx = 0.0000001;
double fEpsGrad = 0.0000001;
int i;
double fBestF;
string sAnswer;
string sErrorMsg = "";
COptimDFP1 oOpt;
MyFxDelegate MyFx = new MyFxDelegate(Fx3);
SayFxDelegate SayFx = new SayFxDelegate(SayFx3);
oOpt = new COptimDFP1();
Console.WriteLine("Davidson-Fletcher-Powell Method (DFP) Optimization");
Console.WriteLine("Finding the minimum of function:");
Console.WriteLine(SayFx());
Console.Write("Use default input values? (Y/N) ");
sAnswer = Console.ReadLine();
if (sAnswer.ToUpper() == "Y")
{
for (i = 0; i < nNumVars; i++)
{
Console.WriteLine("X({0}) = {1}", i + 1, fX[i]);
Console.WriteLine("Tolerance({0}) = {1}", i + 1, fToler[i]);
}
Console.WriteLine("Function tolerance = {0}", fEpsFx);
Console.WriteLine("Maxumum cycles = {0}", nMaxIter);
}
else
{
for (i = 0; i < nNumVars; i++)
{
fX[i] = GetIndexedDblInput("X", i + 1, fX[i]);
fToler[i] = GetIndexedDblInput("Tolerance", i + 1, fToler[i]);
}
fEpsFx = GetDblInput("Function tolerance", fEpsFx);
nMaxIter = GetIntInput("Maxumum cycles", nMaxIter);
}
Console.WriteLine("******** FINAL RESULTS *************");
fBestF = oOpt.CalcOptim(nNumVars, ref fX, ref fParam, ref fToler, fEpsFx, fEpsGrad, nMaxIter, ref nIter, ref sErrorMsg, MyFx);
if (sErrorMsg.Length > 0)
{
Console.WriteLine("** NOTE: {0} ***", sErrorMsg);
}
Console.WriteLine("Optimum at");
for (i = 0; i < nNumVars; i++)
{
Console.WriteLine("X({0}) = {1}", i + 1, fX[i]);
}
Console.WriteLine("Function value = {0}", fBestF);
Console.WriteLine("Number of iterations = {0}", nIter);
Console.WriteLine();
Console.Write("Press Enter to end the program ...");
Console.ReadLine();
}
static public double GetDblInput(string sPrompt, double fDefInput)
{
string sInput;
Console.Write("{0}? ({1}): ", sPrompt, fDefInput);
sInput = Console.ReadLine();
if (sInput.Trim(null).Length > 0)
{
return double.Parse(sInput);
}
else
{
return fDefInput;
}
}
static public int GetIntInput(string sPrompt, int nDefInput)
{
string sInput;
Console.Write("{0}? ({1}): ", sPrompt, nDefInput);
sInput = Console.ReadLine();
if (sInput.Trim(null).Length > 0)
{
return int.Parse(sInput);
}
else
{
return nDefInput;
}
}
static public double GetIndexedDblInput(string sPrompt, int nIndex, double fDefInput)
{
string sInput;
Console.Write("{0}({1})? ({2}): ", sPrompt, nIndex, fDefInput);
sInput = Console.ReadLine();
if (sInput.Trim(null).Length > 0)
{
return double.Parse(sInput);
}
else
{
return fDefInput;
}
}
static public int GetIndexedIntInput(string sPrompt, int nIndex, int nDefInput)
{
string sInput;
Console.Write("{0}({1})? ({2}): ", sPrompt, nIndex, nDefInput);
sInput = Console.ReadLine();
if (sInput.Trim(null).Length > 0)
{
return int.Parse(sInput);
}
else
{
return nDefInput;
}
}
static public string SayFx1()
{
return "F(X) = 10 + (X(1) - 2) ^ 2 + (X(2) + 5) ^ 2";
}
static public double Fx1(int N, ref double[] X, ref double[] fParam)
{
return 10 + Math.Pow(X[0] - 2, 2) + Math.Pow(X[1] + 5, 2);
}
static public string SayFx2()
{
return "F(X) = 100 * (X(1) - X(2) ^ 2) ^ 2 + (X(2) - 1) ^ 2";
}
static public double Fx2(int N, ref double[] X, ref double[] fParam)
{
return Math.Pow(100 * (X[0] - X[1] * X[1]), 2) + Math.Pow((X[1] - 1), 2);
}
static public string SayFx3()
{
return "F(X) = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2";
}
static public double Fx3(int N, ref double[] X, ref double[] fParam)
{
return X[0] - X[1] + 2 * X[0] * X[0] + 2 * X[0] * X[1] + X[1] * X[1];
}
}
}
Notice that the user-defined functions have accompanying helper functions to display the mathematical expression of the function being optimized. For example, function Fx1 has the helper function SayFx1 to list the function optimized in Fx1. Please observe the following rules::
The program uses the following class to optimize the objective function:
using System.Diagnostics;
using System.Data;
using System.Collections;
using Microsoft.VisualBasic;
using System.Collections.Generic;
using System;
namespace Optim_DFP1
{
public delegate double MyFxDelegate(int nNumVars, ref double[] fX, ref double[] fParam);
public delegate string SayFxDelegate();
public class COptimDFP1
{
MyFxDelegate m_MyFx;
protected double MyFxEx(int nNumVars, ref double[] fX, ref double[] fParam, ref double[] fDeltaX, double fLambda)
{
int i;
double[] fXX = new double[nNumVars];
for (i = 0; i < nNumVars; i++)
{
fXX[i] = fX[i] + fLambda * fDeltaX[i];
}
return m_MyFx(nNumVars, ref fXX, ref fParam);
}
protected double FirstDeriv(int nNumVars, ref double[] fX, ref double[] fParam, int nIdxI)
{
double fXt, h, Fp, Fm;
fXt = fX[nIdxI];
h = 0.01 *(1 + Math.Abs(fXt));
fX[nIdxI] = fXt + h;
Fp = m_MyFx(nNumVars, ref fX, ref fParam);
fX[nIdxI] = fXt - h;
Fm = m_MyFx(nNumVars, ref fX, ref fParam);
fX[nIdxI] = fXt;
return (Fp - Fm) / 2 / h;
}
protected void GetFirstDerives(int nNumVars, ref double[] fX, ref double[] fParam, ref double[] fFirstDerivX)
{
int i;
for (i = 0; i < nNumVars; i++)
{
fFirstDerivX[i] = FirstDeriv(nNumVars, ref fX, ref fParam, i);
}
}
protected bool LinSearch_DirectSearch(int nNumVars, ref double[] fX, ref double[] fParam, ref double fLambda, ref double[] fDeltaX, double fInitStep, double fMinStep)
{
double F1, F2;
F1 = MyFxEx(nNumVars, ref fX, ref fParam, ref fDeltaX, fLambda);
do
{
F2 = MyFxEx(nNumVars, ref fX, ref fParam, ref fDeltaX, fLambda + fInitStep);
if (F2 < F1)
{
F1 = F2;
fLambda += fInitStep;
}
else
{
F2 = MyFxEx(nNumVars, ref fX, ref fParam, ref fDeltaX, fLambda - fInitStep);
if (F2 < F1)
{
F1 = F2;
fLambda -= fInitStep;
}
else
{
// reduce search step size
fInitStep /= 10;
}
}
} while (!(fInitStep < fMinStep));
return true;
}
public double CalcOptim(int nNumVars, ref double[] fX, ref double[] fParam, ref double[] fToler, double fEpsFx, double fEpsGrad, int nMaxIter, ref int nIter, ref string sErrorMsg, MyFxDelegate MyFx)
{
int i;
double fLambda;
bool bStop;
int[] Index = new int[nNumVars];
double[] fGrad1 = new double[nNumVars];
double[] fGrad2 = new double[nNumVars];
double[] g = new double[nNumVars];
double[] S = new double[nNumVars];
double[] d = new double[nNumVars];
double[,] Bmat = new double[nNumVars, nNumVars];
double[,] Mmat = new double[nNumVars, nNumVars];
double[,] Nmat = new double[nNumVars, nNumVars];
double[,] MM1 = new double[nNumVars, nNumVars];
double[,] MM2 = new double[nNumVars, nNumVars];
double[,] MM3 = new double[nNumVars, nNumVars];
double[,] MM4 = new double[nNumVars, nNumVars];
double[,] MM5 = new double[nNumVars, nNumVars];
double[] VV1 = new double[nNumVars];
// set error handler
try
{
m_MyFx = MyFx;
// set matrix B as an indentity
MatrixLibVb.MatIdentity(ref Bmat);
nIter = 0;
// calculate initial gradient
GetFirstDerives(nNumVars, ref fX, ref fParam, ref fGrad1);
// start main loop
do
{
nIter++;
if (nIter > nMaxIter)
{
sErrorMsg = "Iter limits";
break;
}
// calcuate vector S() and reverse it's sign
MatrixLibVb.MatMultVect(ref Bmat, ref fGrad1, ref S);
MatrixLibVb.VectMultValInPlace(ref S, - 1);
// normailize vector S
MatrixLibVb.NormalizeVect(ref S);
fLambda = 0.1;
LinSearch_DirectSearch(nNumVars, ref fX, ref fParam, ref fLambda, ref S, 0.1, 0.00001);
// calculate optimum fX() with the given fLambda
for (i = 0; i < nNumVars; i++)
{
d[i] = fLambda * S[i];
fX[i] += d[i];
}
// get new gradient
GetFirstDerives(nNumVars, ref fX, ref fParam, ref fGrad2);
// calculate vector g() and shift fGrad2() into fGrad1()
for (i = 0; i < nNumVars; i++)
{
g[i] = fGrad2[i] - fGrad1[i];
fGrad1[i] = fGrad2[i];
}
// test for convergence
bStop = true;
for (i = 0; i < nNumVars; i++)
{
if (Math.Abs(d[i]) > fToler[i])
{
bStop = false;
break;
}
}
if (bStop)
{
sErrorMsg = "Lam * S()";
break;
}
// exit if gradient is low
if (MatrixLibVb.VectNorm(ref g, 2) < fEpsGrad)
{
sErrorMsg = "G Norm";
break;
}
//-------------------------------------------------
// Start elaborate process to update matrix B
//
// First calculate matrix M (stored as Mmat)
MatrixLibVb.VectToColumnMat(ref S, ref MM1); // MM1 = S in matrix form
MatrixLibVb.VectToRowMat(ref S, ref MM2); // MM2 = S^T in matrix form
MatrixLibVb.VectToColumnMat(ref g, ref MM3); // MM3 = g in matrix form
MatrixLibVb.MatMultMat(ref MM1, ref MM2, ref Mmat); // Mmat = S * S^T
MatrixLibVb.MatMultMat(ref MM2, ref MM3, ref MM4); // MM4 = S^T * g
MatrixLibVb.MatMultValInPlace(ref Mmat, fLambda / MM4[0, 0]); // calculate natrix M
// Calculate matrix nNumVars (stored as Nmat)
MatrixLibVb.MatMultVect(ref Bmat, ref g, ref VV1); // VV1 = {B] g
MatrixLibVb.VectToColumnMat(ref VV1, ref MM1); // MM1 = VV1 in column matrix form ([B] g)
MatrixLibVb.VectToRowMat(ref VV1, ref MM2); // MM2 = VV1 in row matrix form ([B] g)^T
MatrixLibVb.MatMultMat(ref MM1, ref MM2, ref Nmat); // Nmat = ([B] g) * ([B] g)^T
MatrixLibVb.VectToRowMat(ref g, ref MM1); // MM1 = g^T
MatrixLibVb.VectToColumnMat(ref g, ref MM2); // MM2 = g
MatrixLibVb.MatMultMat(ref MM1, ref Bmat, ref MM3); // MM3 = g^T [B]
MatrixLibVb.MatMultMat(ref MM3, ref MM2, ref MM4); // MM4 = g^T [B] g
MatrixLibVb.MatMultValInPlace(ref Nmat, - 1 / MM4[0, 0]); // calculate matrix nNumVars
MatrixLibVb.MatAddMatInPlace(ref Bmat, ref Mmat); // B = B + M
MatrixLibVb.MatAddMatInPlace(ref Bmat, ref Nmat); // B = B + nNumVars
} while (true);
}
catch(Exception ex)
{
sErrorMsg = "Error: " + ex.Message;
}
return MyFx(nNumVars, ref fX, ref fParam);
}
}
}
Copyright (c) Namir Shammas. All rights reserved.