Posted By

mikhailv on 01/31/15


Tagged


Versions (?)

PSO library


 / Published in: C++
 

First experiment in creating a library that enables optimisation by the particle swarm optimisation algorithm

  1. // Class for carrying out optimisation using the particle
  2. // swarm optimization algorithm (PSO)
  3.  
  4. #include<fstream>
  5. #include<string>
  6. #include<iostream>
  7. #include<vector>
  8. #include<cmath>
  9. #include<ctime>
  10. using namespace std;
  11.  
  12.  
  13. // void debug_print(const vector<double> vec); // DEBUGGING
  14.  
  15.  
  16. class udm
  17. { // User-defined methods and constants
  18. public:
  19.  
  20. // Prototype of the objective function
  21. // Precondition: parameter values in the same order as in the input file
  22. // Postcondition: objective function value
  23. static double Objective(const vector<double>& parameterValues);
  24.  
  25. // Prototype of the constraint function
  26. // Precondition: parameter values in the same order as in the input file
  27. // Postcondition: constraints values, (0,1] = infeasible, [-1,0] = feasible
  28. // If no constraints applied the function returns an empty vector
  29. static vector<double> Constraints(const vector<double>& parameterValues);
  30.  
  31. static const int c1 = 4; // local learning factor
  32. static const int c2 = 4; // global learning factor
  33. static const int NOG = 500; // assumed number of generations for calculating inertia weight
  34. };
  35.  
  36.  
  37. class Feasible_set
  38. { // Not to be used by user
  39. public:
  40. Feasible_set()
  41. {
  42. m_numberOfSets = 0;
  43. }
  44.  
  45. // Adds new parameter range at the end of m_data
  46. void add_set(const vector<double>& newset)
  47. {
  48. unsigned int size = newset.size();
  49. for (unsigned int i = 0; i < size; i++)
  50. {
  51. m_data.push_back(newset[i]);
  52. }
  53. m_numbersInEachSet[m_numberOfSets] = size;
  54. m_numberOfSets++;
  55. }
  56.  
  57. // Retrieve parameter value (value number from 1)
  58. double get_value(int setNumber, int valueNumber) const
  59. {
  60. int position = -1; // because m_data starting index is 0
  61. for (int i = 0; i < setNumber - 1; i++)
  62. {
  63. position += m_numbersInEachSet[i];
  64. }
  65. position += valueNumber;
  66. return m_data[position];
  67. }
  68.  
  69. // Retrieve set length (one of the sets)
  70. int get_length(int setnumber) const
  71. {
  72. return m_numbersInEachSet[setnumber - 1];
  73. }
  74.  
  75. // Get the number of dimensions
  76. int numberOfRanges() const
  77. {
  78. return m_numberOfSets;
  79. }
  80.  
  81. private:
  82. vector<double> m_data;
  83. int m_numbersInEachSet[500];
  84. int m_numberOfSets;
  85. };
  86.  
  87.  
  88. class Pso
  89. { // To be used by user
  90. public:
  91. Pso()
  92. { }
  93.  
  94. // Full constructor, takes input file name/path and number of particles
  95. Pso(string filename, unsigned int numberOfPart) :m_inputSource(filename), m_numberOfParticles(numberOfPart)
  96. {
  97. m_seed = 0;
  98. m_readIn(filename);
  99. }
  100.  
  101. // Initialize particles
  102. void initialize()
  103. {
  104. m_initializeAll();
  105. }
  106.  
  107. // Runs new generation for all particles
  108. void operator ++ (int)
  109. {
  110. m_currentGeneration++;
  111.  
  112. for (unsigned int i = 0; i < m_numberOfParticles; i++)
  113. {
  114. m_iterate(m_theParticles[i]);
  115. }
  116. m_determineGBest();
  117. }
  118.  
  119. friend ostream& operator <<(ostream& out, const Pso& opt);
  120.  
  121. // Outputs current global best
  122. double bestFitness() const
  123. {
  124. return m_currentGBest;
  125. }
  126.  
  127. // Outputs current best position
  128. vector<int> bestPosision() const
  129. {
  130. return m_gbestPosition;
  131. }
  132.  
  133. // Accessor to position of a particle
  134. vector<int> currentPosition(unsigned int particleNumber) const
  135. {
  136. return m_theParticles[particleNumber-1].position;
  137. }
  138.  
  139. // Outputs global best parameter values
  140. vector<double> bestValues() const
  141. {
  142. vector<double> values(m_numberOfDimensions);
  143. for (unsigned int i = 0; i < m_numberOfDimensions; i++)
  144. { // Retrieving parameter values at the current position
  145. values[i] = m_myranges.get_value(i + 1, m_gbestPosition[i]);
  146. }
  147. return values;
  148. }
  149.  
  150. private:
  151. struct particle
  152. {
  153. vector<int>position;
  154. vector<int>velocity;
  155. vector<int>pbestPosition;
  156. double pbest;
  157. };
  158.  
  159.  
  160. //////////////// Private variables declarations ////////////////
  161. string m_inputSource;
  162.  
  163. unsigned int m_numberOfDimensions;
  164. unsigned int m_numberOfParticles;
  165.  
  166. int m_currentGeneration;
  167. double m_currentGBest;
  168. vector<int> m_gbestPosition;
  169. vector<int> m_maxVelocityAllowed;
  170.  
  171. vector<particle> m_theParticles;
  172. Feasible_set m_myranges;
  173. int m_seed;
  174.  
  175.  
  176. //////////////// Private functions definitions ////////////////
  177.  
  178. // Reads data ranges into Feasible_set 'm_myranges'
  179. void m_readIn(string filename)
  180. {
  181. ifstream sourceFile;
  182. sourceFile.open(filename, ios::in);
  183. if (!sourceFile.is_open()) {
  184. cerr << "Error: Input file cannot be open\n";
  185. exit(1);
  186. }
  187. m_myranges = Feasible_set();
  188. vector<double> datavec;
  189. int numberOfValues, value;
  190. while (sourceFile >> numberOfValues)
  191. { // reading a line
  192. for (int i = 0; i < numberOfValues; i++)
  193. {
  194. sourceFile >> value;
  195. datavec.push_back(value);
  196. }
  197. m_myranges.add_set(datavec);
  198. datavec.clear();
  199. }
  200. m_numberOfDimensions = m_myranges.numberOfRanges();
  201. }
  202.  
  203. // Calculates inertia weight at the moment
  204. double m_inertiaWeight()
  205. {
  206. return ((udm::NOG*udm::NOG - (m_currentGeneration - udm::NOG)*(m_currentGeneration - udm::NOG))
  207. / (2 * (-udm::NOG*udm::NOG - 2 * udm::NOG + 1)) + 1);
  208. }
  209.  
  210. // Runs one generation for one particle
  211. void m_iterate(particle& mypart)
  212. {
  213. srand(m_seed);
  214. double r1 = rand() / static_cast<double>(RAND_MAX);
  215. double r2 = rand() / static_cast<double>(RAND_MAX);
  216. for (unsigned int i = 0; i < m_numberOfDimensions; i++)
  217. { // Setting new velocity
  218. mypart.velocity[i] = static_cast<int>(m_inertiaWeight()*mypart.velocity[i] + udm::c1*r1*(m_gbestPosition[i] - mypart.position[i])
  219. + udm::c2*r2*(mypart.pbestPosition[i] - mypart.position[i])+0.5);
  220. if (mypart.velocity[i] > m_maxVelocityAllowed[i])
  221. { // Avoids explosions
  222. mypart.velocity[i] = m_maxVelocityAllowed[i];
  223. }
  224. }
  225. for (unsigned int i = 0; i < m_numberOfDimensions; i++)
  226. { // Setting new position
  227. mypart.position[i] += mypart.velocity[i];
  228. while ( mypart.position[i] > m_myranges.get_length(i + 1) )
  229. { // Wall collisions
  230. mypart.position[i]--;
  231. mypart.velocity[i]--;
  232. }
  233. while (mypart.position[i] < 1)
  234. { // Wall collisions
  235. mypart.position[i]++;
  236. mypart.velocity[i]++;
  237. }
  238. }
  239. vector<double> values(m_numberOfDimensions);
  240. for (unsigned int i = 0; i < m_numberOfDimensions; i++)
  241. { // Retrieving parameter values at the current position
  242. values[i] = m_myranges.get_value(i + 1, mypart.position[i]);
  243. }
  244. vector<double> constraints = udm::Constraints(values);
  245.  
  246. double objective = udm::Objective(values);
  247. double fitness = objective;
  248. for (unsigned int i = 0; i < constraints.size(); i++)
  249. { // Final penalty function value
  250. if (constraints[i] > 0)
  251. {
  252. fitness += exp(constraints[i]) * fabs(objective);
  253. }
  254. }
  255. if (fitness < mypart.pbest)
  256. { // New local best
  257. mypart.pbest = fitness;
  258. mypart.pbestPosition = mypart.position;
  259. }
  260. m_seed++;
  261. }
  262.  
  263. void m_determineGBest()
  264. {
  265. for (unsigned int i = 0; i < m_numberOfParticles; i++)
  266. {
  267. if (m_theParticles[i].pbest < m_currentGBest)
  268. { // New global best
  269. m_currentGBest = m_theParticles[i].pbest;
  270. m_gbestPosition = m_theParticles[i].pbestPosition;
  271. }
  272. }
  273. }
  274.  
  275. void m_initializeAll()
  276. {
  277. m_currentGeneration = 1;
  278.  
  279. for (unsigned int i = 1; i <= m_numberOfDimensions; i++)
  280. {
  281. m_maxVelocityAllowed.push_back(static_cast<int>(0.5*m_myranges.get_length(i)));
  282. }
  283.  
  284. vector<particle> temp_part(m_numberOfParticles);
  285. vector<int> temp_dim(m_numberOfDimensions);
  286. for (unsigned int i = 0; i < m_numberOfParticles; i++)
  287. {
  288. temp_part[i].position = temp_dim;
  289. temp_part[i].velocity = temp_dim;
  290. }
  291. m_theParticles = temp_part;
  292.  
  293. double r1;
  294. int r2;
  295. for (unsigned int i = 0; i < m_numberOfParticles; i++)
  296. { // Initial random velocity and position
  297. for (unsigned int j = 0; j < m_numberOfDimensions; j++)
  298. {
  299. srand(m_seed);
  300. r1 = rand() / static_cast<double>(RAND_MAX);
  301. r2 = rand();
  302. m_theParticles[i].velocity[j] = static_cast<int>(m_maxVelocityAllowed[j] * r1 *pow(-1, rand()));
  303. m_theParticles[i].position[j] = r2 % m_myranges.get_length(j + 1) + 1;
  304. m_seed++;
  305. }
  306. }
  307. // Calculating fitness
  308. vector<double> values(m_numberOfDimensions);
  309. for (unsigned int i = 0; i < m_numberOfParticles; i++)
  310. {
  311. for (unsigned int j = 0; j < m_numberOfDimensions; j++)
  312. { // Retrieving parameter values at the current position
  313. values[j] = m_myranges.get_value(j + 1, m_theParticles[i].position[j]);
  314. }
  315.  
  316. // debug_print(values); // DEBUGGING
  317.  
  318. vector<double> constraints = udm::Constraints(values);
  319. double objective = udm::Objective(values);
  320. double fitness = objective;
  321. for (unsigned int i = 0; i < constraints.size(); i++)
  322. { // Final penalty function value
  323. if (constraints[i] > 0)
  324. {
  325. fitness += exp(constraints[i]) * fabs(objective);
  326. }
  327. }
  328. constraints.clear();
  329. m_theParticles[i].pbest = fitness;
  330. m_theParticles[i].pbestPosition = m_theParticles[i].position;
  331. }
  332. // Setting global best
  333. m_currentGBest = m_theParticles[0].pbest;
  334. m_gbestPosition = m_theParticles[0].pbestPosition;
  335. m_determineGBest();
  336. }
  337. };
  338.  
  339. ostream& operator <<(ostream& out, const Pso& opt)
  340. {
  341. out << "Current generation: " << opt.m_currentGeneration << "\n"
  342. << "Current global best: " << opt.m_currentGBest << "\n"
  343. << "Global best parameter values: ";
  344. vector<double> values = opt.bestValues();
  345. for (unsigned int i = 0; i < opt.m_numberOfDimensions; i++)
  346. {
  347. out << values[i] << " ";
  348. }
  349. out << endl << endl;
  350. return out;
  351. }
  352.  
  353.  
  354. // DEBUGGING
  355. //void debug_print(const vector<double> vec)
  356. //{
  357. // cout << endl << endl;
  358. // if (vec.size() == 0)
  359. // cout << "Empty vector";
  360. // else
  361. // {
  362. // for (unsigned int i = 0; i < vec.size(); i++)
  363. // {
  364. // cout << vec[i] << " ";
  365. // }
  366. // }
  367. // cout << endl << endl;
  368. //}

Report this snippet  

You need to login to post a comment.