Hi Dave,
How can I write a routine to minimize the following function using powell method?

Function 1:

eta=0.4;
Doub zRef=-1-1.5/pow(x[0],x[1]);
Doub zHs = (1+eta+SQR(eta)-pow(eta,3))/pow((1-eta),3);


f=(1+x[0]*(zHs-1)+(x[0]-1)*(1+zRef*(x[0]-1)/x[0])*eta*(2*eta-5)/SQR(1-eta)-zFact)

Function 2:

f=10*SQR((SQR(x[1])*(3-x[0])-SQR(x[0])*(3+x[0]))+SQR(2+x[0])/(1+SQR(2+x[0]))

...minimize...using powell method...

#include "../code/nr3.h"
#include "../code/mins.h"
#include "../code/mins_ndim.h"

Doub f(VecDoub_I & x)
{
Doub retval;

//Calculate return value, retval, as a function of x[0], x[1], ...

return retval;
}

int main()
{
VecDoub x(2); // Make its dimension whatever your function is defined for
VecDoub xx; // Its size will be set in the return from minimize()

x[0] = //Initial guess
x[1] = //Initial guess

// Feed the function to the constructor
Powell <Doub(VecDoub_I &)> powell(f);

// Use the minimize member function
xx = powell.minimize(x);

// Print the following if you are interested
//cout << "Number of iterations = " << powell.iter << endl << endl;

cout << scientific;

cout << "Minimum function value of f = "
<< powell.fret
<< ", found at: " << endl;
for (Int i = 0; i < xx.size(); i++) {
cout << " x[" << i << "] = " << setw(13) << xx[i] << endl;
}
cout << endl << endl;

return 0;
}



Regards,

Dave

Thanks Dave,
Here is my code for function 1:

//MixSummary written by Amir Vahid 9/7/2009 in order to do the optimization based on single
//packing fraction
//The optimization code is taken from NR forum 9/22/09
#define _USESTDVECTOR_
#include "C:\Program Files\Numerical Recipes\NR_C302\code\nr3.h"
#include "C:\Program Files\Numerical Recipes\NR_C302\code\mins.h"
#include "C:\Program Files\Numerical Recipes\NR_C302\code\mins_ndim.h"

//Function for calculating of z value based on EG correlation JCP 2005
Doub func2(VecDoub_I & x)
{
Doub MW, Veff;
ifstream DmdWellsSum ("DmdWellsSum.txt");
string txt;
getline(DmdWellsSum,txt);
DmdWellsSum >> MW >> Veff;
ifstream RefTzu ("RefTzu.txt");
VecDoub Eta(32),TempTarget(32),AvgT(32), Tstdev(32), CPUTime(32), density(32),zFactor(32),Zstdev(32), U(32),Udev(32),vEffec(32),vEffstdev(32),Etta(32);
getline(RefTzu,txt);
Int i=0;
Doub eta,zFact;
while (!RefTzu.eof()){
RefTzu >> TempTarget[i]>> AvgT[i] >> Tstdev[i]>> CPUTime[i] >> density [i] >> zFactor[i]>> Zstdev[i] >> U[i] >> Udev[i] >> vEffec[i] >> vEffstdev[i]>> Etta[i];
getline(RefTzu,txt);
Eta[i]=density[i]/MW*(Veff*602.22);
if(abs(Etta[i]-0.2)<0.001){
eta=Eta[i];
zFact=zFactor[i];
}
i=i+1;
}
//Doub zRef;
//zRef=0-1-1.5/pow(x[0],x[1]);
Doub zRef=-1-1.5/pow(x[0],x[1]);
Doub zHs = (1+eta+SQR(eta)-pow(eta,3))/pow((1-eta),3);
return sqrt(pow((1+x[0]*(zHs-1)+(x[0]-1)*(1+zRef*(x[0]-1)/x[0])*eta*(2*eta-5)/SQR(1-eta)-zFact),2));
}
Doub func4(VecDoub_I & y)
{
Doub MW, Veff;
ifstream DmdWellsSum ("DmdWellsSum.txt");
string txt;
getline(DmdWellsSum,txt);
DmdWellsSum >> MW >> Veff;
ifstream RefTzu ("RefTzu.txt");
VecDoub Eta(32),TempTarget(32),AvgT(32), Tstdev(32), CPUTime(32), density(32),zFactor(32),Zstdev(32), U(32),Udev(32),vEffec(32),vEffstdev(32),Etta(32);
getline(RefTzu,txt);
Int i=0;
Doub eta,zFact;
while (!RefTzu.eof()){
RefTzu >> TempTarget[i]>> AvgT[i] >> Tstdev[i]>> CPUTime[i] >> density [i] >> zFactor[i]>> Zstdev[i] >> U[i] >> Udev[i] >> vEffec[i] >> vEffstdev[i]>> Etta[i];
getline(RefTzu,txt);
Eta[i]=density[i]/MW*(Veff*602.22);
if(abs(Etta[i]-0.4)<0.001){
eta=Eta[i];
zFact=zFactor[i];
}
i=i+1;
}
//Doub zRef=0-1-1.5/pow(y[0],y[1]);
Doub zHs = (1+eta+SQR(eta)-pow(eta,3))/pow((1-eta),3);
Doub zRef=-1-1.5/pow(y[0],y[1]);
return sqrt(pow((1+y[0]*(zHs-1)+(y[0]-1)*(1+zRef*(y[0]-1)/y[0])*eta*(2*eta-5)/SQR(1-eta)-zFact),2));
}



int main()
{
Int n = 2; // Number of parameters in the function

Int i;

//
// In general, you might have a better idea for an initial guess, but
// here, I'll just use zero values for the parameters.
//
VecDoub pinit(n, 3); // Initial guess is (3)
VecDoub pmin; // Vector to hold the result from Powell minimization
pinit[1]=-1;
cout << scientific;

cout << "Initial guess: minimum value = " << setw(13) << func2(pinit) << " at: ";
for (i = 0; i < pinit.size(); i++) {
cout << setw(13) << pinit[i] << " ";
}
cout << endl << endl;

VecDoub qinit(n, 3); // Initial guess is (3)
VecDoub qmin; // Vector to hold the result from Powell minimization
qinit[1]=-1;
cout << scientific;
cout << "Initial guess: minimum value = " << setw(13) << func4(qinit) << " at: ";
for (i = 0; i < pinit.size(); i++) {
cout << setw(13) << qinit[i] << " ";
}
cout << endl << endl;

//
// Feed the function to the constructor
//
Powell <Doub(VecDoub_I &)> powell2(func2);

// Perform the minimization
pmin = powell2.minimize(pinit);
cout << "Number of Powell iterations = " << powell2.iter << endl << endl;
cout << "After Powell, minimum value = ";
cout << setw(13) << powell2.fret << " at: ";
for (i = 0; i < n; i++) {
cout << setw(13) << pmin[i] << " ";
}

cout << endl << endl;

Powell <Doub(VecDoub_I &)> powell4(func4);

qmin = powell4.minimize(qinit);
cout << "Number of Powell iterations = " << powell4.iter << endl << endl;
cout << "After Powell, minimum value = ";
cout << setw(13) << powell4.fret << " at: ";
for (i = 0; i < n; i++) {
cout << setw(13) << qmin[i] << " ";
}



cout << endl << endl;

Int nEta=0;
VecDoub TempTarg(32);
ifstream RefTzzu("RefTzu.txt");
string txt;
getline(RefTzzu,txt);
while (!RefTzzu.eof()){
RefTzzu >> TempTarg[nEta];
getline(RefTzzu,txt);
nEta=nEta+1;
}
nEta=nEta-1;


Doub z1,z2,z3, Mw, vEff, zero;
VecDoub eta(32);
ifstream DmdWellsSum ("DmdWellsSum.txt");
//string txt;
getline(DmdWellsSum,txt);
DmdWellsSum >> Mw >> vEff >> zero;
DmdWellsSum >> z1 >> z2 >> z3;
Doub c1,c2,c3;
VecDoub A0(32);
c1 = z1+3;
c2 = z2-3;
c3 = z3+1;
fstream A0SimCalc;
A0SimCalc.open("A0SimCalc.txt", ios::out);
string IP;
VecDoub TempTarget(32),AvgT(32), Tstdev(32), CPUTime(32), density(32), aHs(32),a0Calc2(32),a0Calc4(32);
ifstream RefTzu ("RefTzu.txt");
getline(RefTzu,txt);
cout << "m at eta = 0.2 " << pmin[0] << " " << qmin[0] << endl;
A0SimCalc << "m at eta = 0.2 " << pmin[0] << " " <<pmin[1] << " " << "m at eta = 0.4 " << qmin[0] << " " << qmin[1] << endl;
Doub pRef=-1-1.5/pow(pmin[0],pmin[1]);
Doub qRef=-1-1.5/pow(qmin[0],qmin[1]);
for (int i = 0; i < nEta; i ++){
RefTzu >> TempTarget[i]>> AvgT[i] >> Tstdev[i]>> CPUTime[i] >> density[i];
getline(RefTzu,txt);
eta[i]=density[i]/Mw*(vEff*602.22);
aHs [i] = eta[i]*(4-3*eta[i])/SQR(1-eta[i]);
a0Calc2 [i] = aHs[i]+(pmin[0]-1)*(aHs[i]+(1+pRef*(pmin[0]-1)/pmin[0])*(2*log(1-eta[i])-3*eta[i]/(1-eta[i])));
a0Calc4 [i] = aHs[i]+(qmin[0]-1)*(aHs[i]+(1+qRef*(qmin[0]-1)/qmin[0])*(2*log(1-eta[i])-3*eta[i]/(1-eta[i])));
A0[i] =(c1+c2+c3)/(2*(1-eta[i])*(1-eta[i]))-(c2+2*c3)/(1-eta[i])-c3*log(1-eta[i])-(c1-c2-3*c3)/2;
cout <<eta [i]<< " " << A0[i] << " " << a0Calc2[i] << " " << a0Calc4[i] << endl;
A0SimCalc << eta [i]<< " " << A0[i] << " " << a0Calc2[i] << " " << a0Calc4[i] << endl;
}
A0SimCalc.close();

return 0;
}

I also have enclosed my data files. As you see it only optimizes the first parameter and leaves the 2nd parm at its initial guess. I wonder if you could elaborate on this issue.
Thanks

Here is my code for function 2:

#define _USESTDVECTOR_
#include "C:\Program Files\Numerical Recipes\NR_C302\code\nr3.h"
#include "C:\Program Files\Numerical Recipes\NR_C302\code\mins.h"
#include "C:\Program Files\Numerical Recipes\NR_C302\code\mins_ndim.h"

Doub f(VecDoub_I & x)
{
Doub retval;

retval = 10*SQR((SQR(x[1])*(3-x[0])-SQR(x[0])*(3+x[0]))+SQR(2+x[0])/(1+SQR(2+x[0])));

return retval;
}

int main()
{
VecDoub x(2); // Make its dimension whatever your function is defined for
VecDoub xx; // Its size will be set in the return from minimize()

x[0] = 0.1;//Initial guess
x[1] = 4.2;//Initial guess

// Feed the function to the constructor
Powell <Doub(VecDoub_I &)> powell(f);

// Use the minimize member function
xx = powell.minimize(x);

// Print the following if you are interested
//cout << "Number of iterations = " << powell.iter << endl << endl;

cout << scientific;

cout << "Minimum function value of f = "
<< powell.fret
<< ", found at: " << endl;
for (Int i = 0; i < xx.size(); i++) {
cout << " x[" << i << "] = " << setw(13) << xx[i] << endl;
}
cout << endl << endl;

return 0;
}


As you see this routine does not work either!

...As you see...

I don't see any such thing. As a matter of fact, I don't see anything at all other than code.

Here's the drill:

Tell us what the program is supposed to do. You gave the thread a title that indicates a question about Powell minimization. The program that you posted appears to do more than just Powell stuff.


For purposes of getting help, I suggest that you reduce the program to the part that is problematic. In other words, if your question has to do with Powell stuff, make a smaller program that just does Powell stuff. If your original program does two Powell minimizations, make the test program with only one. I don't think it's reasonable to expect other people to wade through a couple of hundred of lines of other stuff to try to see what the heck you are asking about.


Show us exactly what you expect to see in the output.
Of course you might not know in advance what the exact output values are, but show us what you expect the output to look like. See Footnote.


Show us exactly what you got from the program.


Tell us what exactly you don't understand about the output.


Make the output on your program (the original program and the smaller test program) verbose enough so that you (and we) can know exactly what the program is telling us.

For example, If I might make the output of the test program look something like


Initial guess for func2:
value = 1.913058e+00 at: 3.000000e+00 -1.000000e+00

Number of Powell iterations for func2 = 2

After Powell minimization of func2:
value = 1.084202e-18 at: 1.009404e+00 8.811571e+01


or some such thing.




Regards,

Dave

Footnote: If you don't know (in advance) what the output values are supposed to be, how do you know whether they are correct? How do you test the output to see if its values are reasonable?

Okay. Let's start with function 2 which is easier and has a known minimum at
(x0,x1)=(-2,±0.89442719). This function has been taken from stoer and Bulirsch, Introduction to Numerical Analysis, p.234 and from Numerical Recipes Example Book (C++), p.178. The initial guess is 0.1, 4.2.

I will rewrite the simple code as follows:

#define _USESTDVECTOR_
#include "C:\Program Files\Numerical Recipes\NR_C302\code\nr3.h"
#include "C:\Program Files\Numerical Recipes\NR_C302\code\mins.h"
#include "C:\Program Files\Numerical Recipes\NR_C302\code\mins_ndim.h"

Doub f(VecDoub_I & x)
{
Doub retval;

Doub x1p2sqr=SQR(2.0+x[0]);
retval = 10.0*SQR(SQR(x[1])*(3.0-x[0])-SQR(x[0])*(3.0+x[0]))+x1p2sqr/(1.0+x1p2sqr);

return retval;
}

int main()
{
VecDoub x(2); // Make its dimension whatever your function is defined for
VecDoub xx; // Its size will be set in the return from minimize()

x[0] = 0.1;//Initial guess
x[1] = 4.2;//Initial guess
cout << "Initial guess: minimum value = " << setw(13) << f(x) << " at: ";
for (int i = 0; i < x.size(); i++) {
cout << setw(13) << x[i] << " ";
}
cout << endl << endl;

// Feed the function to the constructor
Powell <Doub(VecDoub_I &)> powell(f);

// Use the minimize member function
xx = powell.minimize(x);

// Print the following if you are interested
//cout << "Number of iterations = " << powell.iter << endl << endl;

cout << scientific;

cout << "Minimum function value of f = "
<< powell.fret
<< ", found at: " << endl;
for (Int i = 0; i < xx.size(); i++) {
cout << " x[" << i << "] = " << setw(13) << xx[i] << endl;
}
cout << endl << endl;

return 0;
}


Here are the initial guess and results:

Initial guess: minimum value = 26138.5 at: 0.1 4.2

Minimum function value of f = 9.395559e-001, found at:
x[0] = 1.942615e+000
x[1] = 4.199992e+000


Press any key to continue . . .


Infact, the results should be something like this:
True minimum is at (-2.0,±0.89442719)
Initial guess: ( 0.1000, 4.2000)

Solution vector: (-2.000000, 0.894427)
Minimum function value 0.000000


I wonder if you could elaborate on this issue. How can the procedure report both +0.8944271 and -0.8944271 in the solution vector?
Warm wishes

I was reading your last reply but it is just deleted!

Okay. Let's start with function 2 which is easier and has a known minimum...

I think that this may be a good example (bad example?) of a function that is not particularly suited for minimization by Powell.

Some methods use gradient information supplied by the user to help them know something about directions to try. It may be that, since Powell is only working with the function and is not given any information about the gradient, it is relatively easier to get to a "wrong" answer with a function that has multiple (and deep) local minima.

Maybe we can get some insight by evaluating the function for a few points close to a solution returned by Powell::minimize.

For example, If we start with initial estimates of 0.1 and 4.2:


Initial guess:
f(x,y) = f(1.00000000000000e-01, 4.20000000000000e+00) = 2.69777992128162e+04

After Powell:
f(x,y) = f(1.96564925319503e+00, 4.19999999484383e+00) = 6.53217159064973e-30

Hold y constant, vary x slightly:
f(x,y) = f(1.96563925319503e+00, 4.19999999484383e+00) = 1.68073664071816e-06
f(x,y) = f(1.96564425319503e+00, 4.19999999484383e+00) = 4.20185073052553e-07
f(x,y) = f(1.96564925319503e+00, 4.19999999484383e+00) = 6.53217159064973e-30
f(x,y) = f(1.96565425319503e+00, 4.19999999484383e+00) = 4.20186898804288e-07
f(x,y) = f(1.96565925319503e+00, 4.19999999484383e+00) = 1.68075124677019e-06

Hold y constant, vary x slightly:
f(x,y) = f(1.96564925319503e+00, 4.19998999484383e+00) = 7.54906564112482e-08
f(x,y) = f(1.96564925319503e+00, 4.19999499484383e+00) = 1.88726865711978e-08
f(x,y) = f(1.96564925319503e+00, 4.19999999484383e+00) = 6.53217159064973e-30
f(x,y) = f(1.96564925319503e+00, 4.20000499484383e+00) = 1.88727315070292e-08
f(x,y) = f(1.96564925319503e+00, 4.20000999484383e+00) = 7.54910158894731e-08


That's a pretty good looking minimum at the values returned by the Powell::minimize routine.

If we try initial values of 0.1 0.894427192587362



Initial guess:
f(x,y) = f(1.00000000000000e-01, 8.94427192587362e-01) = 9.63579145473994e+01

After Powell:
f(x,y) = f(8.28263473076199e-01, 8.94427212930958e-01) = 5.38329392902575e-33

Hold y constant, vary x slightly:
f(x,y) = f(8.28253473076199e-01, 8.94427212930958e-01) = 6.01825116396167e-08
f(x,y) = f(8.28258473076199e-01, 8.94427212930958e-01) = 1.50457348953978e-08
f(x,y) = f(8.28263473076199e-01, 8.94427212930958e-01) = 5.38329392902575e-33
f(x,y) = f(8.28268473076199e-01, 8.94427212930958e-01) = 1.50459488699870e-08
f(x,y) = f(8.28273473076199e-01, 8.94427212930958e-01) = 6.01842234292341e-08

Hold y constant, vary x slightly:
f(x,y) = f(8.28263473076199e-01, 8.94417212930958e-01) = 1.50924385355560e-08
f(x,y) = f(8.28263473076199e-01, 8.94422212930958e-01) = 3.77313072625271e-09
f(x,y) = f(8.28263473076199e-01, 8.94427212930958e-01) = 5.38329392902575e-33
f(x,y) = f(8.28263473076199e-01, 8.94432212930958e-01) = 3.77317291182006e-09
f(x,y) = f(8.28263473076199e-01, 8.94437212930958e-01) = 1.50927760170975e-08


The value here at the point found by Powell is less than the other one. Of course, there is no guarantee (from the Powell stuff) that this is the actual minimum.

In the examples that you cite, you note that Powell goes to a point that has approximately the y value of the initial estimate. With these steep-sided wells, it seems that it could be easy to miss a minimum value while poking around.

Methods that use gradient information from the user, (dfpmin, for example) might be more suitable for certain functions.


How can the procedure report both +0.8944271 and -0.8944271 in the solution vector?
I'm not sure what you are asking here. All of the methods work like the following:



You give it an initial estimate/guess.


The algorithm churns some numbers and returns a point at which it thinks there is a minimum.


You scratch your head and try to figure whether the results are meaningful.



They aren't guaranteed to find an actual minimum, but each one "does its best" and then returns the point that it found.


Maybe if you give it an estimate closer to the other minimum it might settle there (or maybe not).

If you want to see how this implementation looks for minima, just put a print statement in your function definition:

Doub f(VecDoub_I & x)
{
Doub retval;

retval = 10*SQR(
(SQR(x[1])*(3-x[0])-SQR(x[0])*(3+x[0]))+
SQR(2+x[0])/(1+SQR(2+x[0]))
);
cout << "f(" << x[0] << "," << x[1] << ") = " << retval << endl;
return retval;

}



It keeps the y-value (that's x[ 1 ]) constant and varies x (that's x[ 0 ]). When it finds a significant dip, it tries different values of y. Since the steep-sided wells are so narrow, it doesn't find anything worth pursing away from that value of y, so it goes back and changes x and y a little at a time to "zero in" on the sweet spot that it found.


Regards,

Dave

Here is my results. However, I don't know how to get more sig figs from VecDoub and Doub. I wonder if you could elaborate on that. Thanks

f(0.1,4.2) = 26138.5
Initial guess: minimum value = 26138.5 at: 0.1 4.2

f(0.1,4.2) = 26138.5
f(0.1,4.2) = 26138.5
f(1.1,4.2) = 8154.79
f(2.71803,4.2) = 13891
f(1.1,4.2) = 8154.79
f(1.71803,4.2) = 755.735
f(2.1,4.2) = 438.526
f(1.94629,4.2) = 1.16217
f(1.9348,4.2) = 1.94352
f(1.94247,4.2) = 0.939893
f(1.94261,4.2) = 0.939557
f(1.94262,4.2) = 0.939556
f(1.94262,4.2) = 0.939556
f(1.94262,4.2) = 0.939556
f(1.94262,4.2) = 0.939556
f(1.94262,4.2) = 0.939556
f(1.94262,5.2) = 988.866
f(1.94262,2.58197) = 1347.26
f(1.94262,4.2) = 0.939556
f(1.94262,3.58197) = 259.56
f(1.94262,4.58197) = 126.747
f(1.94262,4.17077) = 1.60894
f(1.94262,4.19875) = 0.940787
f(1.94262,4.2) = 0.939556
f(1.94262,4.2) = 0.939556
f(1.94262,4.2) = 0.939556
f(1.94262,4.2) = 0.939556
f(1.94262,4.2) = 0.939556
f(1.94262,4.2) = 0.939556
f(1.94262,4.2) = 0.939556
f(1.94262,4.2) = 0.939556
f(1.94262,4.2) = 0.939556
f(3.78523,4.19999) = 123367
f(1.94262,4.2) = 0.939556
f(2.94262,4.2) = 25447.7
f(0.324582,4.2) = 21944.5
f(1.94262,4.2) = 0.939556
f(1.32458,4.2) = 4826.31
f(2.32458,4.2) = 2842.86
f(1.88962,4.2) = 46.2137
f(1.93852,4.2) = 1.21585
f(1.94266,4.2) = 0.939594
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,4.2) = 0.939556
f(1.94261,5.2) = 988.866
f(1.94261,2.58196) = 1347.26
f(1.94261,4.2) = 0.939556
f(1.94261,3.58196) = 259.56
f(1.94261,4.58196) = 126.747
f(1.94261,4.17076) = 1.60894
f(1.94261,4.19874) = 0.940787
f(1.94261,4.2) = 0.939556
f(1.94261,4.19999) = 0.939556
f(1.94261,4.19999) = 0.939556
f(1.94261,4.19999) = 0.939556
f(1.94261,4.19999) = 0.939556
f(1.94261,4.19999) = 0.939556
f(1.94261,4.19999) = 0.939556
f(1.94261,4.19999) = 0.939556
Minimum function value of f = 9.395559e-001, found at:
x[0] = 1.942615e+000
x[1] = 4.199992e+000


Press any key to continue . . .

...more sig figs from...Doub....

I often use scientific notation for floating point numbers since that's the simplest way to keep things lined up (and for other even more important reasons: namely the treatment of numbers with small absolute value is uniform). Use setprecision() to tell it how many digits to print after the decimal point:


cout << scientific << setprecision(15);
cout << "22/700 = " << 22.0/700.0 << endl;




Output

22/700 = 3.142857142857143e-02


These parameters hold for the rest of the program unless you want to change them. In other words, you don't need to do them before each and every output statement.

Regards,

Dave