Begin implementation of generic PSO optimization tools

CMakeLists.txt

@@ -20,6 +20,7 @@ include_directories(./src/Quadrature/include/)
 include_directories(./src/Elements/include/)
 include_directories(./src/Matrix/include/)
 include_directories(./src/Materials/include/)
+include_directories(./src/Optimization/include/)
 include_directories(./src/Physics/include/)
 include_directories(./src/Model\ Templates/include/)
 

src/CMakeLists.txt

@@ -5,6 +5,7 @@ include_directories(./Quadrature/)
 include_directories(./Elements/)
 include_directories(./Matrix/)
 include_directories(./Materials/)
+include_directories(./Optimization/)
 include_directories(./Physics/)
 include_directories(./Model\ Templates/)
 
@@ -15,5 +16,6 @@ add_subdirectory(./Quadrature/)
 add_subdirectory(./Elements/)
 add_subdirectory(./Matrix/)
 add_subdirectory(./Materials/)
+add_subdirectory(./Optimization/)
 add_subdirectory(./Physics/)
 add_subdirectory(./Model\ Templates/)
\ No newline at end of file

src/Optimization/CMakeLists.txt

+project(Oersted_Optimization)
+
+set(SOURCE_FILES
+        ./include/Optimization.hpp
+
+        include/Optimization.hpp src/Particle.cpp src/Particle.h)
+
+add_library(optimization SHARED ${SOURCE_FILES})
+
+target_link_libraries(optimization)
\ No newline at end of file

src/Optimization/include/Optimization.hpp

+#ifndef OERSTED_OPTIMIZATION_HPP
+#define OERSTED_OPTIMIZATION_HPP
+
+#include "../src/Particle.h"
+
+#endif //OERSTED_OPTIMIZATION_HPP

src/Optimization/src/Particle.cpp

+#include "Particle.h"
+
+void Particle::update_personal_best(std::map<std::string,double_t> go) {
+        double_t value = LocalObjective(go);
+
+        if (value < PersonalBest) {
+            HasImproved = true;
+            PersonalBest = value;
+            for (auto &key_state : State) {
+                auto const& key = key_state.first;
+                personal_best(key) = position(key);
+            }
+        }
+}
+
+void Particle::update_global_best(std::map<std::string,double_t> go, Particle &p) {
+    double_t value = LocalObjective(go);
+
+    if (value < GlobalBest) {
+        HasImproved = true;
+        GlobalBest = value;
+        for (auto &key_state : State) {
+            std::string const &key = key_state.first;
+            global_best(key) = p.position(key);
+        }
+    }
+}
+
+void Particle::update_velocity(std::mt19937 rng, double_t omega, double_t phip, double_t phig) {
+    std::uniform_real_distribution<> dist(0.0,1.0);
+
+    for (auto &key_state : State) {
+        CoordinateState &state = key_state.second;
+        state.update_velocity(omega, dist(rng) * phip, dist(rng) * phig);
+    }
+}
+
+void Particle::update_position() {
+    for (auto &key_state : State) {
+        CoordinateState &state = key_state.second;
+        state.update_position();
+    }
+}
+
+bool Particle::has_improved() {
+    bool output = HasImproved;
+    HasImproved = false;
+    return output;
+}
+
+bool Particle::is_converged(CoordinateSpace &space, double_t state_tol, double_t objective_tol) {
+    if (PersonalBest - GlobalBest > objective_tol) {
+        return false;
+    }
+
+    for (auto &key_state : State) {
+        std::string key = key_state.first;
+        CoordinateState &state = key_state.second;
+
+        if(!state.is_converged(state_tol * space.coordinate(key).range())) {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+std::ostream& operator<<(std::ostream& os, const Particle& p) {
+    os << "P-Best: " << p.PersonalBest << ", G-Best: " << p.GlobalBest;
+}
+
+Particle CoordinateSpace::new_particle(std::mt19937 &rng) {
+    std::uniform_real_distribution<> dist(-1.0,1.0);
+
+    Particle particle;
+    for (auto &key_coord : Coordinates) {
+        std::string const &key = key_coord.first;
+        Coordinate const &coord = key_coord.second;
+
+        double_t v = dist(rng) * coord.range() * SpeedLimit;
+        double_t p = dist(rng) * coord.range() * SpeedLimit + coord.nominal();
+
+        particle.state(key) = CoordinateState{p,v};
+    }
+
+    return particle;
+}
+
+void CoordinateSpace::limit_velocity(Particle &p) {
+    for (auto &key_coord : Coordinates) {
+        std::string const &key = key_coord.first;
+        Coordinate const &coord = key_coord.second;
+
+        double_t &pv = p.velocity(key);
+        double_t v = SpeedLimit * coord.range();
+        if (pv > v) {
+            pv = v;
+        } else if (pv < -v) {
+            pv = -v;
+        }
+    }
+}
+
+void CoordinateSpace::limit_position(Particle &p) {
+    for (auto &key_coord : Coordinates) {
+        std::string const &key = key_coord.first;
+        Coordinate const &coord = key_coord.second;
+
+        double_t &px = p.position(key);
+        double_t &pv = p.velocity(key);
+
+        double_t lb = coord.lower_bound();
+        double_t ub = coord.upper_bound();
+        if (px < lb) {
+            pv = (lb - px);
+            px = lb;
+        } else if (px > ub) {
+            pv = (ub - px);
+            px = ub;
+        }
+    }
+}
+
+void CoordinateSpace::perturb(Particle &p, std::mt19937 &rng){
+    std::uniform_real_distribution<> dist(-Perturbation,+Perturbation);
+    for (auto &key_coord : Coordinates) {
+        std::string const &key = key_coord.first;
+        Coordinate const &coord = key_coord.second;
+
+        p.position(key) += dist(rng) * coord.range();
+        p.velocity(key) += dist(rng) * coord.range();
+    }
+}
+
+void Swarm::run() {
+    std::uniform_real_distribution<> dist(0.0,1.0);
+
+    bool all_converged{false};
+    for (size_t iter = 0; iter != MaximumIterations && !all_converged; ++iter) {
+        std::cout << iter << std::endl;
+        for (size_t i = 0; i != Particles.size(); ++i) {
+            auto go = GlobalObjective(Particles[i]);
+
+            Particles[i].update_personal_best(go);
+            for (size_t j = 0; j != Particles.size(); ++j) {
+                Particles[j].update_global_best(go, Particles[i]);
+            }
+        }
+
+        all_converged = true;
+        for (auto &p : Particles) {
+            all_converged &= p.is_converged(Space,CoordinateTolerance,ObjectiveTolerance);
+        }
+
+        if (!all_converged) {
+            for (auto &p : Particles) {
+                if (!p.has_improved()) {
+                    Space.perturb(p,RNG);
+                }
+            }
+
+            Omega = 2.0 / (1 + M_SQRT2) * dist(RNG);
+
+            for (auto &p : Particles) {
+                p.update_velocity(RNG, Omega, PhiP, PhiG);
+
+                Space.limit_velocity(p);
+
+                p.update_position();
+
+                Space.limit_position(p);
+            }
+        }
+    }
+}
\ No newline at end of file

src/Optimization/src/Particle.h

+#ifndef OERSTED_PARTICLE_H
+#define OERSTED_PARTICLE_H
+
+#include <cmath>
+#include <cfloat>
+#include <map>
+#include <string>
+#include <functional>
+#include <random>
+#include <iostream>
+
+class Particle;
+class CoordinateSpace;
+class Coordinate;
+class CoordinateState;
+
+class CoordinateState {
+public:
+    CoordinateState() {};
+
+    CoordinateState(double_t p, double_t v) : Position{p}, Velocity{v}, PersonalBest{p}, GlobalBest{p} {};
+
+    void update_velocity(double_t omega, double_t phip, double_t phig) { Velocity = omega * (Velocity + phip * (PersonalBest - Position) + phig * (GlobalBest - Position)); }
+
+    void update_position() { Position += Velocity; };
+
+    bool is_converged(double_t tol) { return std::abs(PersonalBest - GlobalBest) < tol; };
+
+    double_t Position;
+
+    double_t Velocity;
+
+    double_t PersonalBest;
+
+    double_t GlobalBest;
+};
+
+class Coordinate {
+public:
+    Coordinate() {};
+
+    Coordinate(double_t lb, double_t nom, double_t ub) : LowerBound{lb}, Nominal{nom}, UpperBound{ub} {};
+
+    double_t lower_bound() const { return LowerBound; };
+
+    double_t nominal() const { return Nominal; };
+
+    double_t upper_bound() const { return UpperBound; };
+
+    double_t range() const { return UpperBound - LowerBound; };
+
+protected:
+    double_t LowerBound;
+    double_t Nominal;
+    double_t UpperBound;
+};
+
+class Particle {
+public:
+    CoordinateState &state(std::string key) { return State[key]; };
+
+    double_t &position(std::string key) { return State[key].Position; };
+
+    double_t &velocity(std::string key) { return State[key].Velocity; };
+
+    double_t &personal_best(std::string key) { return State[key].PersonalBest; };
+
+    double_t &global_best(std::string key) { return State[key].GlobalBest; };
+
+    double_t &personal_best() { return PersonalBest; };
+
+    double_t &global_best() { return GlobalBest; };
+
+    std::function<double_t(std::map<std::string,double_t>)> &local_objective() { return LocalObjective; };
+
+    double_t local_objective(std::map<std::string,double_t> &go) { return LocalObjective(go); };
+
+    void update_personal_best(std::map<std::string,double_t> objectives);
+
+    void update_global_best(std::map<std::string,double_t> objectives, Particle &p);
+
+    void update_velocity(std::mt19937 rng, double_t omega, double_t phip, double_t phig);
+
+    void update_position();
+
+    bool has_improved();
+
+    bool is_converged(CoordinateSpace &ss, double_t coordinate_tol, double_t objective_tol);
+
+    friend std::ostream& operator<<(std::ostream& os, const Particle& p);
+
+    std::function<double_t(std::map<std::string,double_t>)> LocalObjective;
+
+    double_t PersonalBest{DBL_MAX};
+
+    double_t GlobalBest{DBL_MAX};
+
+protected:
+    std::map<std::string,CoordinateState> State;
+
+    bool HasImproved{true};
+
+};
+
+
+class CoordinateSpace {
+public:
+    Coordinate &coordinate(std::string key) { return Coordinates[key]; };
+
+    Particle new_particle(std::mt19937 &rng);
+
+    size_t size() const { return Coordinates.size(); };
+
+    void limit_velocity(Particle &p);
+
+    void limit_position(Particle &p);
+
+    void perturb(Particle &p, std::mt19937 &rng);
+
+    std::map<std::string,Coordinate> Coordinates;
+
+    double_t SpeedLimit{0.2};
+
+    double_t Perturbation{0.2e-2};
+};
+
+class Swarm {
+public:
+    Swarm(CoordinateSpace ss, std::function<std::map<std::string,double_t>(Particle)> go, size_t Np, double_t otol, size_t maxit) : Space{ss}, GlobalObjective{go}, ObjectiveTolerance{otol}, MaximumIterations{maxit} {
+        std::random_device rd;
+        RNG = std::mt19937{rd()};
+
+        for (size_t i = 0; i != Np; ++i) {
+            Particles.push_back(Space.new_particle(RNG));
+        }
+    }
+
+    size_t size() const { return Particles.size(); };
+
+    void run();
+
+    std::vector<Particle> Particles;
+
+    double_t Omega{1.0 / (1.0 + M_SQRT2)};
+    double_t PhiP{4.0};
+    double_t PhiG{4.0};
+
+    double_t CoordinateTolerance{1e-2};
+    double_t ObjectiveTolerance;
+    size_t MaximumIterations;
+
+protected:
+    CoordinateSpace Space;
+
+    std::function<std::map<std::string,double_t>(Particle)> GlobalObjective;
+
+    std::mt19937 RNG;
+};
+
+#endif //OERSTED_PARTICLE_H

test/CMakeLists.txt

@@ -30,6 +30,9 @@ set(SOURCE_FILES
         Physics/Finite_Element_Mesh_Test.cpp
         Physics/test_SquareOrder2.cpp
 
+        Optimization/test_Optimization.hpp
+        Optimization/test_ParticleSwarmOptimization.cpp
+
         LibraryIntegration/Library_Integration_Test.hpp
         LibraryIntegration/Mesh_To_FEM_Test.cpp
         LibraryIntegration/test_Square.cpp
@@ -43,8 +46,8 @@ set(SOURCE_FILES
 
         Model\ Templates/test_ModelTemplates.hpp
         Model\ Templates/test_DistributedWindingStator.cpp
-        Model\ Templates/test_SynchronousReluctanceRotor.cpp "Model Templates/test_ModelTemplates.hpp" "Model Templates/test_SynchronousReluctanceMachine.cpp")
+        Model\ Templates/test_SynchronousReluctanceRotor.cpp "Model Templates/test_ModelTemplates.hpp" "Model Templates/test_SynchronousReluctanceMachine.cpp" Optimization/test_Optimization.hpp Optimization/test_ParticleSwarmOptimization.cpp)
 
 add_executable(run_tests ${SOURCE_FILES})
 
-target_link_libraries(run_tests gtest gtest_main geometry sketch elements mesh materials quadrature physics model_templates stdc++fs)
\ No newline at end of file
+target_link_libraries(run_tests gtest gtest_main geometry sketch elements mesh materials optimization quadrature physics model_templates stdc++fs)
\ No newline at end of file

test/Optimization/test_Optimization.hpp

+#ifndef OERSTED_TEST_OPTIMIZATION_HPP
+#define OERSTED_TEST_OPTIMIZATION_HPP
+
+#include "Optimization.hpp"
+#include "gtest.h"
+
+#endif //OERSTED_TEST_OPTIMIZATION_HPP

test/Optimization/test_ParticleSwarmOptimization.cpp

+#include "test_Optimization.hpp"
+
+TEST(Particle_Swarm_Optimization, coordinate_space_initialization) {
+    CoordinateSpace cs;
+
+    cs.coordinate("x") = Coordinate{-1.0,0.0,1.0};
+    cs.coordinate("y") = Coordinate{2.0,3.0,5.0};
+
+    EXPECT_EQ(cs.coordinate("x").lower_bound(), -1.0);
+    EXPECT_EQ(cs.coordinate("x").nominal(), 0.0);
+    EXPECT_EQ(cs.coordinate("x").upper_bound(), +1.0);
+    EXPECT_NEAR(cs.coordinate("x").range(), 2.0, FLT_EPSILON);
+
+
+    EXPECT_EQ(cs.coordinate("y").lower_bound(), 2.0);
+    EXPECT_EQ(cs.coordinate("y").nominal(), 3.0);
+    EXPECT_EQ(cs.coordinate("y").upper_bound(), 5.0);
+    EXPECT_NEAR(cs.coordinate("y").range(), 3.0, FLT_EPSILON);
+}
+
+TEST(Particle_Swarm_Optimization, particle_initialization) {
+    CoordinateSpace cs;
+
+    cs.coordinate("x") = Coordinate{-1.0,0.0,1.0};
+    cs.coordinate("y") = Coordinate{2.0,3.0,5.0};
+    
+    std::random_device rd;
+    std::mt19937 rng(rd());
+
+    Particle p = cs.new_particle(rng);
+
+    for (std::string key : {"x","y"}) {
+        EXPECT_GE(p.position(key), cs.coordinate(key).lower_bound());
+        EXPECT_LE(p.position(key), cs.coordinate(key).upper_bound());
+
+        EXPECT_GE(p.velocity(key), -cs.coordinate(key).range());
+        EXPECT_LE(p.velocity(key), +cs.coordinate(key).range());
+
+        EXPECT_GE(p.personal_best(key), cs.coordinate(key).lower_bound());
+        EXPECT_LE(p.personal_best(key), cs.coordinate(key).upper_bound());
+
+        EXPECT_GE(p.global_best(key), cs.coordinate(key).lower_bound());
+        EXPECT_LE(p.global_best(key), cs.coordinate(key).upper_bound());
+
+        EXPECT_GE(p.personal_best(), FLT_MAX);
+        EXPECT_GE(p.global_best(), FLT_MAX);
+    }
+
+    double_t w0{1.0/3.0};
+    double_t w1{2.0/3.0};
+    p.local_objective() = [w0,w1](std::map<std::string,double_t> go) {
+        return w0*go["z0"]+w1*go["z1"];
+    };
+
+    std::map<std::string,double_t> go;
+    go["z0"] = M_PI;
+    go["z1"] = M_E;
+
+    EXPECT_NEAR(p.local_objective(go), w0 * go["z0"] + w1 * go["z1"], FLT_EPSILON);
+}
+
+TEST(Particle_Swarm_Optimization, single_objective_swarm) {
+    CoordinateSpace cs;
+
+    cs.coordinate("x") = Coordinate{-2.0,0.0,2.0};
+    cs.coordinate("y") = Coordinate{-2.0,0.0,2.0};
+
+    std::function<std::map<std::string,double_t>(Particle p)> go = [](Particle p) {
+        std::map<std::string,double_t> output;
+
+        double_t x = p.position("x");
+        double_t y = p.position("y");
+
+        output["z0"] = (x - 0.5) * (x - 0.5) + (y - 0.5) * (y - 0.5);
+
+        return output;
+    };
+
+    Swarm swarm{cs,go,10,1e-2,20};
+
+    for (size_t i = 0; i != swarm.size(); ++i) {
+        swarm.Particles[i].local_objective() = [](std::map<std::string,double_t> go) {
+            return go["z0"];
+        };
+    }
+
+    swarm.run();
+
+    for (size_t i = 0; i != swarm.size(); ++i) {
+        std::cout << i << ": " << swarm.Particles[i] << std::endl;
+        std::cout << swarm.Particles[i].personal_best("x") << "," << swarm.Particles[i].personal_best("y") << std::endl;
+    }
+
+    GTEST_FATAL_FAILURE_("TODO: IMPLEMENT TESTS");
+    GTEST_FATAL_FAILURE_("TODO: CLEAN FILE STRUCTURE");
+}
+
+TEST(Particle_Swarm_Optimization, multiobjective_swarm) {
+    CoordinateSpace cs;
+
+    cs.coordinate("x") = Coordinate{-2.0,0.0,2.0};
+    cs.coordinate("y") = Coordinate{-2.0,0.0,2.0};
+
+    std::function<std::map<std::string,double_t>(Particle p)> go = [](Particle p) {
+        std::map<std::string,double_t> output;
+
+        double_t x = p.position("x");
+        double_t y = p.position("y");
+
+        output["z0"] = (x - 0.5) * (x - 0.5) + (y - 0.5) * (y - 0.5);
+        output["z1"] = (x + 0.5) * (x + 0.5) + (y - 0.5) * (y - 0.5);
+
+        return output;
+    };
+
+    Swarm swarm{cs,go,11,0.25e-2,30};
+
+    for (size_t i = 0; i != swarm.size(); ++i) {
+        double_t w0 = (10 - i) / 10.0;
+        double_t w1 = i / 10.0;
+
+        swarm.Particles[i].local_objective() = [w0,w1](std::map<std::string,double_t> go) {
+            return w0 * go["z0"] + w1 * go["z1"];
+        };
+    }
+
+    swarm.run();
+
+    for (size_t i = 0; i != swarm.size(); ++i) {
+        std::cout << i << ": " << swarm.Particles[i] << std::endl;
+        std::cout << swarm.Particles[i].personal_best("x") << "," << swarm.Particles[i].personal_best("y") << std::endl;
+    }
+
+    GTEST_FATAL_FAILURE_("TODO: IMPLEMENT TESTS");
+    GTEST_FATAL_FAILURE_("TODO: CLEAN FILE STRUCTURE");
+}
\ No newline at end of file