The following program calculates the minimum point of a multi-variable function using the Powell search method.
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 for each variable:
1. The function tolerance value
2. The step tolerance value.
3. The maximum number of search cycles
4. The values for the initial set of variables
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 value.
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 an initial value of 0 for each variable, a function tolerance of 1e-7, and a maximum number of 100 search cycles. 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_PowellSearch1
{
sealed class Module1
{
static public void Main()
{
int nNumVars = 2;
double[] fX = new double[] { 0, 0 };
double[] fParam = new double[] { 0, 0 };
int nIter = 0;
int nMaxCycles = 100;
double fEpsFx = 0.0000001;
double fEpsStep = 0.0000001;
int i;
double fBestF;
string sAnswer;
string sErrorMsg = "";
CPowellSearch1 oOpt;
MyFxDelegate MyFx = new MyFxDelegate( Fx3);
SayFxDelegate SayFx = new SayFxDelegate( SayFx3);
oOpt = new CPowellSearch1();
Console.WriteLine("Powell Search Optimization (version 1)");
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("Function tolerance = {0}", fEpsFx);
Console.WriteLine("Step tolerance = {0}", fEpsStep);
Console.WriteLine("Maxumum cycles = {0}", nMaxCycles);
}
else
{
for (i = 0; i < nNumVars; i++)
{
fX[i] = GetIndexedDblInput("X", i + 1, fX[i]);
}
fEpsFx = GetDblInput("Function tolerance", fEpsFx);
fEpsStep = GetDblInput("Function tolerance", fEpsStep);
nMaxCycles = (int) GetDblInput("Maxumum cycles", nMaxCycles);
}
Console.WriteLine("******** FINAL RESULTS *************");
fBestF = oOpt.CalcOptim(nNumVars, ref fX, ref fParam, fEpsFx, fEpsStep, nMaxCycles, ref nIter, ref sErrorMsg, MyFx);
if (sErrorMsg.Length > 0)
{
Console.WriteLine("**ERROR: {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) double.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) double.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_PowellSearch1
{
public delegate double MyFxDelegate(int nNumVars, ref double[] fX, ref double[] fParam);
public delegate string SayFxDelegate();
public class CPowellSearch1
{
MyFxDelegate m_MyFx;
protected double MyFxEx(int nNumVars, ref double[] fX, ref double[] fParam, ref double[,] fSmat, int nII, double fAlpha)
{
int i;
double[] fXX = new double[nNumVars + 1];
for (i = 0; i < nNumVars; i++)
{
fXX[i] = fX[i] + fAlpha * fSmat[nII, i];
}
return m_MyFx(nNumVars, ref fXX, ref fParam);
}
protected bool LinSearch_DirectSearch(int nNumVars, ref double[] fX, ref double[] fParam, ref double fAlpha, ref double[,] fSmat, int nII, double fInitStep, double fMinStep)
{
double F1, F2;
F1 = MyFxEx(nNumVars, ref fX, ref fParam, ref fSmat, nII, fAlpha);
do
{
F2 = MyFxEx(nNumVars, ref fX, ref fParam, ref fSmat, nII, fAlpha + fInitStep);
if (F2 < F1)
{
F1 = F2;
fAlpha += fInitStep;
}
else
{
F2 = MyFxEx(nNumVars, ref fX, ref fX, ref fSmat, nII, fAlpha - fInitStep);
if (F2 < F1)
{
F1 = F2;
fAlpha -= fInitStep;
}
else
{
// reduce search step size
fInitStep /= 10;
}
}
} while (!(fInitStep < fMinStep));
return true;
}
public double CalcOptim(int nNumVars, ref double[] fX, ref double[] fParam, double fEpsFx, double fEpsStep, int nMaxCycles, ref int nIter, ref string sErrorMsg, MyFxDelegate MyFx)
{
int i, j, k, m;
double fAlpha, fSum, F1, F2;
double[] fX1 = new double[nNumVars];
double[] fX2 = new double[nNumVars];
double[,] fSmat = new double[nNumVars + 1, nNumVars];
m_MyFx = MyFx;
F1 = MyFx(nNumVars, ref fX, ref fParam);
// Build fSmat matrix
for (k = 0; k < nNumVars; k++)
{
fX1[k] = fX[k];
for (m = 0; m < nNumVars; m++)
{
fSmat[k, m] = k == m ? 1 : 0;
}
}
sErrorMsg = "";
nIter = 1;
j = 1;
do
{
nIter++;
// reset row fSmat(nNumVars+1,:)
for (k = 0; k < nNumVars; k++)
{
fSmat[nNumVars, k] = 0;
}
for (i = 0; i < nNumVars; i++)
{
fAlpha = 0;
if (!LinSearch_DirectSearch(nNumVars, ref fX, ref fParam, ref fAlpha, ref fSmat, i, 0.1, 0.0001))
{
sErrorMsg = "Linear Search failed";
break;
}
for (k = 0; k < nNumVars; k++)
{
fX[k] += fAlpha * fSmat[i, k];
fSmat[nNumVars, k] = fSmat[nNumVars, k] + fAlpha * fSmat[i, k];
}
}
fAlpha = 0;
if (!LinSearch_DirectSearch(nNumVars, ref fX, ref fX, ref fAlpha, ref fSmat, nNumVars, 0.1, 0.0001))
{
sErrorMsg = "Linear Search failed";
break;
}
for (k = 0; k < nNumVars; k++)
{
fX[k] += fAlpha * fSmat[nNumVars, k];
fX2[k] = fX[k];
}
object temp_object4 = fParam;
double[] temp_double4 = (double[]) temp_object4;
F2 = MyFx(nNumVars, ref fX2, ref temp_double4);
// test end of iterations criteria
if (Math.Abs(F2 - F1) < fEpsFx)
{
break;
}
fSum = 0;
for (k = 0; k < nNumVars; k++)
{
fSum += Math.Pow((fX2[k] - fX1[k]), 2);
}
fSum = Math.Sqrt(fSum);
if (fSum < fEpsStep)
{
break;
}
if (j >= nMaxCycles)
{
break;
}
j++;
// rotate data
for (k = 0; k < nNumVars; k++)
{
fX1[k] = fX2[k];
}
F1 = F2;
// copy fSmat matrix
for (k = 0; k < nNumVars; k++)
{
for (m = 0; m < nNumVars; m++)
{
fSmat[k, m] = fSmat[k + 1, m];
}
}
} while (true);
return F1;
}
}
}
Copyright (c) Namir Shammas. All rights reserved.