Tmp
Putting it all together -- Optimization and EPAnet
This section contains some advanced stuff. Don't be discouraged if you don't understand all and difficult to complete it without the help of somebody else.
Prerequisites
- You should have completed the C++ programming primer preceeding this section.
- You have exposure to the EPAnet 2 software (graphical version is adequate).
- A basic understanding of Genetic Algorithms (Follow these links: [1],[2],[3] and spend some time if not.)
In this lesson, we push the abilities we have gained to the limit! We link two code libraries, namely, Evolving objects -- a versatile genetic algorithm code and EPAnet (toolkit) a pipe network calculation model by writing some relatively simple code and create a running program. The following section explains the problem we have selected to solve (please note that the problem itself is of secondary importance here, what we are trying to do is to hone the skills we have gained and build confidence to attack bigger coding projects).
You can download the set of files needed for this exercise from : EO.rar.
Problem
The figure on the left shows a water supply network with a reservoir, seven pipe segments conveying water and five junctions that have different demands. We want to compute the most economical pipe diameters for each segment while maintaining a minimum pressure head of 10 m at all junctions. You can open the File:Cbcpp pipe network1.inp in Epanet 2 software to view the network.
Our plan is to use a Genetic Algorithm to optimize the pipe diameters. For this we need to connect Epanet Software to a genetic algorithm code.
Plan
We will use the Epanet Toolkit, a programming interface designed to run Epanet 2.0 programmatically without the standard graphical interface. For the Genetic Algorithm part, we'll use Evolving objects a free and open source package for evolutionary computations.
Evolving objects can be downloaded from eodev] website. Epanet toolkit can be downloaded from [US-EPA]. However, for this exercise you may use the already downloaded versions of these tools. Download the file cbcpp_EO.zip into your computer and extract it to a new folder (say E:\GA\EO). This will create four subfolders,
code data eo-1.0 epanet_tools
We use code folder to keep all the code we write (not the code libraries we use!). It already have the files:
RealEA.cpp real_value.h
RealEA.cpp is a c++ program that can be used to access the EO Genetic algorithm functions. real_value.h has the cost function for the GA. We'll learn more about this later.
data folder where we shall keep all our data, has the water distribution network file.
eo-1.0 and epanet_tools have the evolving objects and epanet toolkit code libraries.
Solution
We shall attack our problem in a piecemeal fashion, in the steps given below:
- Get EO to solve a small GA problem.
- Replace the cost function with our own.
- Use epanet_toolkit to run Epanet on our water supply network and use the results to evaluate a cost function.
Usually this step-by-step approach is less error-prone than attempting the whole task once.
Running a GA
Create a new folder called projects under the same folder that has sub-folders of code, data, etc. This is where we shall keep the visual studio express edition related (platform specific) stuff. Open visual C++ and reate a new empty project EPGA in that folder. Add the following files in the code folder to the project.
RealEA.cpp real_value.h
When you try to compile the project, you will get the error message on the line:
#include <eo>
Indicating that the compiler can not include 'eo'. To remedy this, we should include the path to the include file eo ("..\..\eo-1.0.1\src"). This can be done by additing it to: Edit-><project>Properties->C/C++->General->Additional include directories.
At this stage RealEA.cpp should comile successfully, but would cause errors at linking stage. The error messages would look like
unresolved external symbol "public: virtual __thiscall eoFunctorStore::~eoFunctorStore(void)" (??1eoFunctorStore@@UAE@XZ) referenced in function "public: virtual void * __thiscall eoFunctorStore::`scalar deleting destructor'(unsigned int)" (??_GeoFunctorStore@@UAEPAXI@Z)
The reason for this is
- The compiler knows about eo library (by #include <eo>), but
- the real library objects of eo needed for linking are missing
To rectify this, we should let the linker access to the necessary libraries. Before doing this we have to make a detour. First save your project/solution and close it.
Building EO libraries
The eo code available for download does not have the libraries for windows built in, but they have the complete source and tools needed to build them. There is a visual C++ solution called eo.sln in the sub-folder win. Just build this solution (Project->Build solution). This will create the necessary libraries in the win.
One more point: It is beneficial to build the release version (not the debug version) of the libraries for performance reasons.
After this you will have the following libraries in the folder win\lib\release.
eo.lib eoes.lib eoga.lib eoutils.lib
Now open your project again. In add the above four files as dependancies for the linker. (Edit-><project>Properties->Linker->General->Additional Dependancies).
Then let the linker know where these are: (Edit-><project>Properties->Linker->General->Additional Library Directories). Add something like ..\..\eo-1.0.1\win\lib\release.
At this stage, you should be able to compile the project successfully. Debug->Start without Debugging should run the program, albeit with not-so-meaningful-at-the-moment results.
Enable debugging
If you try to debug your code at this stage, visual c++ 2005 will complain (as of 01:07, 9 April 2007 (JST)) that the binaries were not built with debug information. There is a separate article on how to enable debugging in visual studio 2005.
EO In-Action
Now is a good time to have an idea about how our GA code works. Don't worry if you can not understand everything -- what is important is to have a general idea of how things work!
The actual code-in-action is very short, indeed. I have changed the comments a little bit.
/* Many instructions to this program can be given on the command line.
The following code understands what you have specified as command line arguments.
*/
eoParser parser(argc, argv); // for user-parameter reading
eoState state; // keeps all things allocated
typedef eoReal<eoMinimizingFitness> EOT;
// The evaluation fn - encapsulated into an eval counter for output
eoEvalFuncPtr<EOT, double, const std::vector<double>&>
mainEval( real_value );
eoEvalFuncCounter<EOT> eval(mainEval);
// the genotype - through a genotype initializer
eoRealInitBounded<EOT>& init = make_genotype(parser, state, EOT());
// Build the variation operator (any seq/prop construct)
eoGenOp<EOT>& op = make_op(parser, state, init);
// initialize the population - and evaluate
// yes, this is representation indepedent once you have an eoInit
eoPop<EOT>& pop = make_pop(parser, state, init);
// stopping criteria
eoContinue<EOT> & term = make_continue(parser, state, eval);
// output
eoCheckPoint<EOT> & checkpoint = make_checkpoint(parser, state, eval, term);
// algorithm (need the operator!)
eoAlgo<EOT>& ea = make_algo_scalar(parser, state, eval, checkpoint, op);
make_help(parser); //print help, if something is missing or user gives /h
// evaluate intial population AFTER help and status in case it takes time
apply<EOT>(eval, pop);
// print it out
cout << "Initial Population\n";
pop.sortedPrintOn(cout);
cout << endl;
run_ea(ea, pop); // run the ea
cout << "Final Population\n";
pop.sortedPrintOn(cout);
cout << endl;
You can get away without understanding a single line of code above!! The only critical part is the following function, specified in the header file real_value.h. In order to adopt these versatile algorithms to solve a problem of our choosing, only changing this header file is adequate.
Fitness function
double real_value(const std::vector<double>& _ind)
{
double sum = 0;
for (unsigned i = 0; i < _ind.size(); i++)
sum += _ind[i] * _ind[i];
return sqrt(sum);
}
The eo library expects this function to be present and uses to evaluate the individuals in the population for fitness.
What happens here?
- Our genotype vector has (say) n number of items. The fitness function simply squares each of these and add them up and returns the square root. No rocket science here! Just, the rule here is larger the individuals -- better the fitness!!
- The GA code does the rest,
Try running the algorithm. It will go on decreasing the cost function and exit at 100th generation. The default is 100 generations, in a moment, we'll get into the details on how to change the default behavior.
Providing arguments
Once you run the above code once successfully, check the folder that consists the executable file (project\EPGA\Debug folder). There should be a file EPGA.exe.status EO library creates this file, when you run the program. It lists all the changes you can make at the runtime to the GA. (Meaning, you don't have to recompile the program -- just indicate the parameters.) There are two ways of specifying parameters -- first is to list them in the command line. For example try the following in the dos prompt.
EPGA.exe --help
It will just print a long help message. Then try
EPGA.exe --maxGen=1000
It will run a GA like before, but will stop only after 1000 generations.
EPGA.exe --maxGen=0 --steadyGen=20
will run the GA until 20 generations without any improvement.
So, you get the idea.
Now try this. First copy the file EPGA.exe.status to args.txt and edit the latter with notepad.
> copy EPGA.exe.status args.txt
> notepad args.txt
Notice that this file has all the arguments that can be given at the command line. # in a line means the rest that follows is a comment. Note that almost all the argments are commented. Now make sure only the following items are uncommented.
--vecSize=3
--maxGen=0
--steadyGen=20
Now, run the program with the following command.
EPGA.exe @args.txt
Notice that this is equivalent to running
EPGA.exe --vecSize=3 --maxGen=0 --steadyGen=20
Within Visual Studio
To make things a bit easier for us, lets do the following. Move the file args.txt to the data folder.
> move args.txt ..\..\data\.
Then in your project in visual C++, add the file data.txt to resource files tab, in solutions explorer. (This step is optional)
Finally, add the following
@..\..\data\args.txt
to the Command argments Configuration properties->Debugging.
Now within visual C++ itself, you can just modify the file args.txt to suit your needs, save it and re-run our program.
Writing our own cost function
Now let's turn into our original problem. We need to minimize the diameter of seven pipes while maintaining reasonable (>10m) pressure at all supply nodes. Lets forget the pressure part for a moment and concentrate on minimizing the diameter. (The answer here is obvious, if pressure is no concern, then the 'best' diameter is zero!! -- but let's live with this silly notion for the moment.)
Do the following changes:
- First we change the args.txt file to have
--vecSize=7 # we have seven pipes, each diameter can be represented by a single real value.
--initBounds=7[0,500] # -B : Bounds for variables -- pipe diameter within 0-500mm
--objectBounds=7[0,500] # -B : Bounds for variables -- pipe diameter within 0-500mm
- Change the cost function real_value.h to the following:
#include <vector>
#include <iostream>
using namespace std;
double real_value(const std::vector<double>& _ind)
{
//GA returns diameter as ind_
double length=1000;
double factor=1.; //some factor so that cost=length=diameter*factor (lets keep things simple!)
double dia,cost = 0;
for (unsigned i = 0; i < _ind.size(); i++){
cost+=_ind[i]*length*factor;
}
return cost/10000;
}
If everything is all right, you will have some large initial cost, and a very small final cost.
Now that we have customized the GA to represent our silly optimization problem, it is but a small step to do the real job!
In comes EPANet 2
With the provided EPANET toolkit (in epanet_tools folder) there is a help file: TOOLKIT.HLP. This is going to be our standard reference to various calls to epanet program via toolkit.
Let's do most of the changes in the real_value.h file and try to keep changes in RealEA.cpp to a minimum.
We will focus on several toolkit functions:
ENOpen and ENClose
- Declaration
int ENopen( char* f1, char* f2, char* f3)
- Description
- Opens the Toolkit to analyze a particular distribution system.
- Declaration
int Enclose(void)
- Description
- Closes down the Toolkit system (including all files being processed).
ENSolveH
- Declaration
int ENsolveH( void )
- Description
- Runs a complete hydraulic simulation with results for all time periods written to the binary Hydraulics file.
ENGetNodeValue
- Declaration
int ENgetnodevalue( int index, int paramcode, float* value )
- Description
- Retrieves the value of a specific link parameter.
ENSetLinkValue
- Declaration
int ENsetlinkvalue( int index, int paramcode, float value )
- Description
- Sets the value of a parameter for a specific link.
Call EPANet 2
Do the following changes in RealEA.cpp.
try
{
// first initialize the Epanet
epanet_init(); // <-- new addition A function call to initialize epanet. we have to write this function.
eoParser parser(argc, argv); // for user-parameter reading
...
// close Epanet
epanet_close(); // <-- new addition A function call to close epanet. we have to write this function.
}
catch(exception& e)
Now let's do the major modifications in real_value.h
As the first stage:
#include <vector>
#include <iostream>
#include "epanet2.h"
using namespace std;
double dia_cost_factor=1.; //some factor so that cost=length=diameter*factor (lets keep things simple!)
/** A function that computes the cost. This is what the GA use to evaluate its populations */
double real_value(const std::vector<double>& _ind)
{
//GA returns diameter as ind_
double length=1000; /* All pipe lengths are equal */
double dia,cost = 0;
for (unsigned i = 0; i < _ind.size(); i++){
cost+=_ind[i]*length*dia_cost_factor;
}
return cost/10000;
}
/* We open the epanet system with the necessary input file.
A bit of hard coding here. But, lets live with that for the moment. */
void epanet_init(){
int ret;
char file[500]="../../data/network.inp";
char rept[500]="../../data/network.rep";
ret=ENopen(file,rept,"");
cout << "At opening Epanet retured : "<<ret<<'\n';
}
/* Close the epanet system */
void epanet_close(){
int ret;
ret=ENclose();
cout << "At closing Epanet retured : "<<ret<<'\n';
}
To run the above you will have to
- Add ..\..\epanet_tools to additional include directories.
- Add epanet2vc.lib to additional dependencies.
- Add ..\..\epanet_tools to additional library directories (At this stage the application should compile, but will give a runtime error.)
- Add PATH=%PATH%;..\..\epanet_tools to debugging environment. (Project properties->Debugging->Environment)
Run the application at this stage to make sure that the two epanet calls return zero (error free call signal).
Then add a function pressure_cost to real_value.h to compute the 'cost' of pressure deficiency. (Something like the one below)
/* Returns the pressure cost (penalty for pressure violations at demand nodes) based on epanet runs.
Prerequisites: The epanet system should be initialized before calling this function for the first time. */
double pressure_cost(vector<double> _ind){
int ret;
double cost;
for(unsigned int i=0;i<npipes;i++){
int index=-1;
ret=ENgetlinkindex(pipes[i],&index);
//cout << "At opening Epanet retured : "<<ret<<'\n';
ret=ENsetlinkvalue(index,EN_DIAMETER,_ind[i]);
//cout << "At opening Epanet retured : "<<ret<<'\n';
}
//now run the simulation
ret=ENsolveH();
//cout << "At solve Epanet retured : "<<ret<<'\n';
cost=0;
//read the pressure values
for(unsigned int i=0;i<nnodes;i++){
int index=-1;
ret=ENgetnodeindex(nodes[i],&index);
float value;
//cout << "At ENgetnodeindex Epanet retured : "<<ret<<'\n';
ret=ENgetnodevalue(index,EN_PRESSURE,&value);
//cout << "At ENgetnodevalue Epanet retured : "<<ret<<'\n';
if(value<10){
cost+=pressue_cost_factor*(10-value); // if p<10m, set a proportional penalty.
}
}
//cout << "At ENcloseH Epanet retured : "<<ret<<'\n';
return cost;
}
The value of variable pressure_cost_factor should be carefully considered (against that of dia_cost_factor).
Finally modify real_value function so that it will call the above pressure_cost function and add the cost to the total cost.
At this stage you have a complete living-breathing program that join the power to evolving objects with Epanet.
A touch of sophistication -- let's get rid of hard coding
This section is here mostly for the sake of completion. Ignore this section if you don't have time or inclination.
We can change the behavior of the GA without recompiling the code, thanks to the sophistication of the design of EO library. However, we have given up some of this flexibility in the way we have designed our cost function. We have hard coded a number of items:
- Name of the network file.
- Number of pipes in the network.
- Number of nodes.
- IDs of pipes and nodes.
It is possible to make our program quite flexible in these aspects also. But it needs a bit of work. Let's see how it can be done. The first stage is to change the real_value.h so that instead of hard coded values, it can take values stored in variables. Then in the main function (RealEA.cpp) add the code necessary to read these from a text file supplied by user.
Lets define the text file format as follows:
<network file name> <file name for the report to write into> <no of pipes> <PIPE_ID1> ... <no of nodes> <NODE_ID1> ...
That means something like the file below
..\..\data\network.inp ..\..\data\network.rpt 7 P1 P2 P3 P4 P5 P6 P7 5 J1 J2 J3 J4 J5
Then the modified program is as follows:
- real_value.h
#include <vector>
#include <iostream>
#include "epanet2.h"
#include <fstream>
using namespace std;
#define MAX_PATH_LEN 500
#define MAX_LABEL_LEN 25
double pressure_cost(vector<double> _ind);
int npipes=7;
int nnodes=5;
char file[MAX_PATH_LEN];
char rept[MAX_PATH_LEN];
vector<string> nodes; // notice, now we use vectors instead of traditional arrays.
vector<string> pipes;
double dia_cost_factor=1.; //some factor so that cost=length=diameter*factor (lets keep things simple!)
double pressue_cost_factor=1000000; //multiplier to 'map' pressue defficiency to cost.
// cost=(pressuredefficiency)*pmult
/** read the text file specified by filename argument and obtain epanet related parameters */
void parse_epanet_para(char* filename){
cout << "I read epanet related data from "<<filename<<"\n"; // inform the user
//open the file
ifstream myfile (filename);
if(!myfile.is_open()){ // this is important.
cout << "I can not open the file:"<<filename <<" I quit!!\n";
exit(1);
}
myfile >> file; //read the name of the file
myfile >> rept; //read the name of the (new) report file
myfile >> npipes; //number of pipes
for(int i=0;i<npipes;i++){ // read those pipe ids
char tmp[MAX_LABEL_LEN];
myfile >> tmp;
pipes.push_back(tmp);
}
myfile >> nnodes; //number of junctions
for(int i=0;i<nnodes;i++){//those ids
char tmp[MAX_LABEL_LEN];
myfile >> tmp;
nodes.push_back(tmp);
}
}
double real_value(const std::vector<double>& _ind)
{
// check for sanity
if(_ind.size()!=npipes){
//raise hell
cout << "Bloody murder!\n";
cout << "Number of pipes and chromosome size mismatch!\n";
exit(5);
}
//GA returns diameter as ind_
double length=1000;
double dia,cost = 0;
cost=pressure_cost(_ind);
for (unsigned i = 0; i < _ind.size(); i++){
cost+=_ind[i]*length*dia_cost_factor;
}
return cost/10000;
}
double pressure_cost(vector<double> _ind){
int ret;
double cost;
for(unsigned int i=0;i<npipes;i++){
int index=-1;
char tmp[MAX_LABEL_LEN]; // this gimmick here is to convet a c++ string to a c style char*
strcpy(tmp,pipes[i].c_str()); // because epanet is writtin in old c, which does not accept strings.
ret=ENgetlinkindex(tmp,&index);
//cout << "At opening Epanet retured : "<<ret<<'\n';
ret=ENsetlinkvalue(index,EN_DIAMETER,_ind[i]);
//cout << "At opening Epanet retured : "<<ret<<'\n';
}
//now run the simulation
ret=ENsolveH();
//cout << "At solve Epanet retured : "<<ret<<'\n';
cost=0;
//read the pressure values
for(unsigned int i=0;i<nnodes;i++){
int index=-1;
char tmp[MAX_LABEL_LEN]; // convert c++ string to c style char*
strcpy(tmp,nodes[i].c_str());
ret=ENgetnodeindex(tmp,&index);
float value;
//cout << "At ENgetnodeindex Epanet retured : "<<ret<<'\n';
ret=ENgetnodevalue(index,EN_PRESSURE,&value);
//cout << "At ENgetnodevalue Epanet retured : "<<ret<<'\n';
if(value<10){
cost+=pressue_cost_factor*(10-value);
}
}
//cout << "At ENcloseH Epanet retured : "<<ret<<'\n';
return cost;
}
void epanet_init(){
int ret;
ret=ENopen(file,rept,"");
cout << "At opening Epanet retured : "<<ret<<'\n';
}
void epanet_close(){
int ret;
ret=ENclose();
cout << "At closing Epanet retured : "<<ret<<'\n';
}
- RealEA.cpp
#define _CRT_SECURE_NO_DEPRECATE
//above is to get rid of deprecation warnings of Microsoft compiler. Needed because we use strcpy() function.
#include <iostream>
#include <es/make_real.h>
#include "real_value.h"
#include <apply.h>
using namespace std;
typedef eoReal<eoMinimizingFitness> EOT;
void print_values(eoPop<EOT> pop);
int main_function(int argc, char* argv[]);
/** Notice that we have moved everything that was previously in main() to
main_function.
Now before GA related stuff is handled (by main_function),
We process the command argument list. Unlike the previous case, now the first argument, i.e.
the filename of the EPAnet related parameters, is mandatory.
Then we copy the rest of the arguments in argv to a new array argv_ and pass it to main_function.
From the GA viewpoint, nothing has changed. It receives a argument array. If there are no arguments in it, GA will
run with default parameters. Otherwise it will parse the argument array.
The first command line argument is separately passed parse_epanet_para function.
*/
int main(int argc,char *argv[]){
/* argv[0] is always the name of the program. So, to run properly the program should have
length of argv (i.e. argc) >=2
If this is not the case, provide some help. */
if(argc<2){// no arguments provided at the command line
cout << "Usage: argv[0] <epanet_related_datafile> <EO related arguments ...>\n";
cout << "Format of epanet_related_datafile as follows.\n";
cout << "<network file name>\n";
cout << "<file name for the report to write into>\n";
cout << "<no of pipes>\n";
cout << "<PIPE_ID1>\n";
cout << "...\n";
cout << "<no of nodes>\n";
cout << "<NODE_ID1>\n";
cout << "...\n";
exit(1);
}
char* filename=argv[1]; // seperately copy argv[1] (first argument) to variable filename
char* argv_[MAX_PATH_LEN];
argv_[0]=argv[0]; // argv[0] is the calling program name, copy this as is.
for(int i=1;i<argc-1;i++){ // then copy the rest (argv[2], argv[3], ...) of arguments to new array
cerr << argv[i+1];
argv_[i]=argv[i+1];
}
argc--; // argc should be one less than before
//now parse the parameter file stright away!
parse_epanet_para(filename);
return main_function(argc,argv_); // now call main_function with new argc, argv_ pair.
}
int main_function(int argc, char* argv[])
{
try
{
// first initialize the Epanet
epanet_init();
eoParser parser(argc, argv); // for user-parameter reading
eoState state;
eoEvalFuncPtr<EOT, double, const std::vector<double>&>
mainEval( real_value );
eoEvalFuncCounter<EOT> eval(mainEval);
eoRealInitBounded<EOT>& init = make_genotype(parser, state, EOT());
// Build the variation operator (any seq/prop construct)
eoGenOp<EOT>& op = make_op(parser, state, init);
// initialize the population - and evaluate
// yes, this is representation indepedent once you have an eoInit
eoPop<EOT>& pop = make_pop(parser, state, init);
// stopping criteria
eoContinue<EOT> & term = make_continue(parser, state, eval);
// output
eoCheckPoint<EOT> & checkpoint = make_checkpoint(parser, state, eval, term);
// algorithm (need the operator!)
eoAlgo<EOT>& ea = make_algo_scalar(parser, state, eval, checkpoint, op);
// to be called AFTER all parameters have been read!!!
make_help(parser);
// evaluate intial population AFTER help and status in case it takes time
apply<EOT>(eval, pop);
// print it out
cout << "Initial Population\n";
pop.sortedPrintOn(cout);
cout << endl;
cin.get();
run_ea(ea, pop); // run the ea
cout << "Final Population\n";
pop.sortedPrintOn(cout);
cout << endl;
// close Epanet
epanet_close();
}
catch(exception& e)
{
cout << e.what() << endl;
}
return 1;
}
