416 lines
17 KiB
C++
416 lines
17 KiB
C++
#include "reducedmodeloptimizer.hpp"
|
|
#include "bobyqa.h"
|
|
#include "gradientDescent.h"
|
|
#include "matplot/matplot.h"
|
|
#include "simulationhistoryplotter.hpp"
|
|
|
|
const bool draw = true;
|
|
size_t numberOfOptimizationRounds{0};
|
|
Simulator simulator;
|
|
Eigen::MatrixX3d optimalReducedModelDisplacements;
|
|
SimulationJob reducedModelSimulationJob;
|
|
std::unordered_map<ReducedModelVertexIndex, FullModelVertexIndex>
|
|
reducedToFullVertexIndexMap;
|
|
matplot::line_handle plotHandle;
|
|
std::vector<double> fValueHistory;
|
|
Eigen::Vector4d initialParameters;
|
|
struct OptimizationCallback {
|
|
double operator()(const size_t &iterations, const Eigen::VectorXd &x,
|
|
const double &fval, Eigen::VectorXd &gradient) const {
|
|
// run simulation
|
|
// SimulationResults reducedModelResults =
|
|
// simulator.executeSimulation(reducedModelSimulationJob);
|
|
// reducedModelResults.draw(reducedModelSimulationJob);
|
|
fValueHistory.push_back(fval);
|
|
auto xPlot =
|
|
matplot::linspace(0, fValueHistory.size(), fValueHistory.size());
|
|
plotHandle = matplot::scatter(xPlot, fValueHistory);
|
|
// const std::string plotImageFilename = "objectivePlot.png";
|
|
// matplot::save(plotImageFilename);
|
|
// if (numberOfOptimizationRounds % 30 == 0) {
|
|
// std::filesystem::copy_file(
|
|
// std::filesystem::path(plotImageFilename),
|
|
// std::filesystem::path("objectivePlot_copy.png"));
|
|
// }
|
|
// std::stringstream ss;
|
|
// ss << x;
|
|
// reducedModelResults.simulationLabel = ss.str();
|
|
// SimulationResultsReporter resultsReporter;
|
|
// resultsReporter.reportResults(
|
|
// {reducedModelResults},
|
|
// std::filesystem::current_path().append("Results"));
|
|
|
|
return true;
|
|
}
|
|
};
|
|
|
|
struct Objective {
|
|
double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const {
|
|
assert(x.rows() == 4);
|
|
|
|
// drawSimulationJob(simulationJob);
|
|
// Set mesh from x
|
|
std::shared_ptr<ElementalMesh> reducedModel =
|
|
reducedModelSimulationJob.mesh;
|
|
for (EdgeIndex ei = 0; ei < reducedModel->EN(); ei++) {
|
|
BeamFormFinder::Element &e = reducedModel->elements[ei];
|
|
e.properties.A = initialParameters(0) * x(0);
|
|
e.properties.J = initialParameters(1) * x(1);
|
|
e.properties.I2 = initialParameters(2) * x(2);
|
|
e.properties.I3 = initialParameters(3) * x(3);
|
|
e.axialConstFactor = e.properties.E * e.properties.A / e.initialLength;
|
|
e.torsionConstFactor = e.properties.G * e.properties.J / e.initialLength;
|
|
e.firstBendingConstFactor =
|
|
2 * e.properties.E * e.properties.I2 / e.initialLength;
|
|
e.secondBendingConstFactor =
|
|
2 * e.properties.E * e.properties.I3 / e.initialLength;
|
|
}
|
|
// run simulation
|
|
SimulationResults reducedModelResults =
|
|
simulator.executeSimulation(reducedModelSimulationJob);
|
|
// std::stringstream ss;
|
|
// ss << x;
|
|
// reducedModelResults.simulationLabel = ss.str();
|
|
// SimulationResultsReporter resultsReporter;
|
|
// resultsReporter.reportResults(
|
|
// {reducedModelResults},
|
|
// std::filesystem::current_path().append("Results"));
|
|
// compute error and return it
|
|
double error = 0;
|
|
for (const auto reducedFullViPair : reducedToFullVertexIndexMap) {
|
|
VertexIndex reducedModelVi = reducedFullViPair.first;
|
|
Eigen::Vector3d vertexDisplacement(
|
|
reducedModelResults.displacements[reducedModelVi][0],
|
|
reducedModelResults.displacements[reducedModelVi][1],
|
|
reducedModelResults.displacements[reducedModelVi][2]);
|
|
Eigen::Vector3d errorVector =
|
|
Eigen::Vector3d(
|
|
optimalReducedModelDisplacements.row(reducedModelVi)) -
|
|
vertexDisplacement;
|
|
error += errorVector.norm();
|
|
}
|
|
|
|
return error;
|
|
}
|
|
};
|
|
|
|
double objective(long n, const double *x) {
|
|
Eigen::VectorXd eigenX(n, 1);
|
|
for (size_t xi = 0; xi < n; xi++) {
|
|
eigenX(xi) = x[xi];
|
|
}
|
|
// Eigen::VectorXd emptyGradient;
|
|
// return operator()(eigenX, emptyGradient);
|
|
|
|
// drawSimulationJob(simulationJob);
|
|
// Set mesh from x
|
|
std::shared_ptr<ElementalMesh> reducedModel = reducedModelSimulationJob.mesh;
|
|
for (EdgeIndex ei = 0; ei < reducedModel->EN(); ei++) {
|
|
BeamFormFinder::Element &e = reducedModel->elements[ei];
|
|
e.properties.A = initialParameters(0) * eigenX(0);
|
|
e.properties.J = initialParameters(1) * eigenX(1);
|
|
e.properties.I2 = initialParameters(2) * eigenX(2);
|
|
e.properties.I3 = initialParameters(3) * eigenX(3);
|
|
e.axialConstFactor = e.properties.E * e.properties.A / e.initialLength;
|
|
e.torsionConstFactor = e.properties.G * e.properties.J / e.initialLength;
|
|
e.firstBendingConstFactor =
|
|
2 * e.properties.E * e.properties.I2 / e.initialLength;
|
|
e.secondBendingConstFactor =
|
|
2 * e.properties.E * e.properties.I3 / e.initialLength;
|
|
}
|
|
// run simulation
|
|
SimulationResults reducedModelResults =
|
|
simulator.executeSimulation(reducedModelSimulationJob);
|
|
// std::stringstream ss;
|
|
// ss << x;
|
|
// reducedModelResults.simulationLabel = ss.str();
|
|
// SimulationResultsReporter resultsReporter;
|
|
// resultsReporter.reportResults(
|
|
// {reducedModelResults},
|
|
// std::filesystem::current_path().append("Results"));
|
|
// compute error and return it
|
|
double error = 0;
|
|
for (const auto reducedFullViPair : reducedToFullVertexIndexMap) {
|
|
VertexIndex reducedModelVi = reducedFullViPair.first;
|
|
Eigen::Vector3d vertexDisplacement(
|
|
reducedModelResults.displacements[reducedModelVi][0],
|
|
reducedModelResults.displacements[reducedModelVi][1],
|
|
reducedModelResults.displacements[reducedModelVi][2]);
|
|
Eigen::Vector3d errorVector =
|
|
Eigen::Vector3d(optimalReducedModelDisplacements.row(reducedModelVi)) -
|
|
vertexDisplacement;
|
|
error += errorVector.norm();
|
|
}
|
|
|
|
return error;
|
|
}
|
|
|
|
ReducedModelOptimizer::ReducedModelOptimizer(
|
|
ConstVCGEdgeMesh &fullModel, ConstVCGEdgeMesh &reducedModel,
|
|
const std::unordered_map<ReducedModelVertexIndex, FullModelVertexIndex>
|
|
&viMap,
|
|
const std::vector<std::pair<ReducedModelVertexIndex,
|
|
ReducedModelVertexIndex>> &oppositeVertices)
|
|
: reducedToFullViMap(viMap) {
|
|
reducedToFullVertexIndexMap = viMap;
|
|
|
|
// construct fullToReducedViMap
|
|
for (const std::pair<ReducedModelVertexIndex, FullModelVertexIndex> &viPair :
|
|
viMap) {
|
|
fullToReducedViMap[viPair.second] = viPair.first;
|
|
}
|
|
|
|
// construct opposite vertex map
|
|
for (const std::pair<VertexIndex, VertexIndex> &viPair : oppositeVertices) {
|
|
oppositeVertexMap[viPair.first] = viPair.second;
|
|
oppositeVertexMap[viPair.second] = viPair.first;
|
|
}
|
|
|
|
reducedModelElementalMesh = std::make_shared<ElementalMesh>(reducedModel);
|
|
pFullModelElementalMesh = std::make_shared<ElementalMesh>(fullModel);
|
|
optimalReducedModelDisplacements.resize(reducedModelElementalMesh->VN(), 3);
|
|
}
|
|
|
|
void ReducedModelOptimizer::computeReducedModelSimulationJob(
|
|
const SimulationJob &simulationJobOfFullModel,
|
|
SimulationJob &simulationJobOfReducedModel) {
|
|
std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
|
|
reducedModelFixedVertices;
|
|
for (auto fullModelFixedVertex : simulationJobOfFullModel.fixedVertices) {
|
|
reducedModelFixedVertices[fullToReducedViMap.at(
|
|
fullModelFixedVertex.first)] = fullModelFixedVertex.second;
|
|
}
|
|
|
|
std::unordered_map<VertexIndex, Vector6d> reducedModelNodalForces;
|
|
for (auto fullModelNodalForce :
|
|
simulationJobOfFullModel.nodalExternalForces) {
|
|
reducedModelNodalForces[fullToReducedViMap.at(fullModelNodalForce.first)] =
|
|
fullModelNodalForce.second;
|
|
}
|
|
|
|
simulationJobOfReducedModel = SimulationJob{reducedModelElementalMesh,
|
|
reducedModelFixedVertices,
|
|
reducedModelNodalForces,
|
|
{}};
|
|
}
|
|
|
|
SimulationJob ReducedModelOptimizer::getReducedSimulationJob(
|
|
const SimulationJob &fullModelSimulationJob) {
|
|
SimulationJob reducedModelSimulationJob;
|
|
computeReducedModelSimulationJob(fullModelSimulationJob,
|
|
reducedModelSimulationJob);
|
|
return reducedModelSimulationJob;
|
|
}
|
|
|
|
void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
|
|
const SimulationResults &fullModelResults,
|
|
Eigen::MatrixX3d &optimaldisplacementsOfReducedModel) {
|
|
for (auto reducedFullViPair : reducedToFullVertexIndexMap) {
|
|
const VertexIndex fullModelVi = reducedFullViPair.second;
|
|
Vector6d fullModelViDisplacements =
|
|
fullModelResults.displacements[fullModelVi];
|
|
optimaldisplacementsOfReducedModel.row(reducedFullViPair.first) =
|
|
Eigen::Vector3d(fullModelViDisplacements[0],
|
|
fullModelViDisplacements[1],
|
|
fullModelViDisplacements[2]);
|
|
}
|
|
}
|
|
|
|
Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob(
|
|
const SimulationJob &fullModelSimulationJob) {
|
|
SimulationResults fullModelResults =
|
|
simulator.executeSimulation(fullModelSimulationJob, false, false);
|
|
fullModelResults.simulationLabel = "fullModel";
|
|
if (draw) {
|
|
drawWorldAxes();
|
|
fullModelResults.draw(fullModelSimulationJob);
|
|
}
|
|
computeDesiredReducedModelDisplacements(fullModelResults,
|
|
optimalReducedModelDisplacements);
|
|
computeReducedModelSimulationJob(fullModelSimulationJob,
|
|
reducedModelSimulationJob);
|
|
|
|
// Set initial guess of solution
|
|
Eigen::VectorXd initialGuess(4);
|
|
auto propertiesOfFirstBeamOfFullModel =
|
|
pFullModelElementalMesh->elements[0].properties;
|
|
initialParameters(0) = propertiesOfFirstBeamOfFullModel.A;
|
|
initialParameters(1) = propertiesOfFirstBeamOfFullModel.J;
|
|
initialParameters(2) = propertiesOfFirstBeamOfFullModel.I2;
|
|
initialParameters(3) = propertiesOfFirstBeamOfFullModel.I3;
|
|
const double stifnessFactor = 2;
|
|
initialGuess(0) = stifnessFactor;
|
|
initialGuess(1) = stifnessFactor;
|
|
initialGuess(2) = stifnessFactor;
|
|
initialGuess(3) = stifnessFactor;
|
|
|
|
const bool useGradientDescent = false;
|
|
if (useGradientDescent) {
|
|
|
|
// gdc::GradientDescent<double, Objective,
|
|
// gdc::DecreaseBacktracking<double>,
|
|
// OptimizationCallback>
|
|
gdc::GradientDescent<double, Objective, gdc::BarzilaiBorwein<double>,
|
|
OptimizationCallback>
|
|
// gdc::GradientDescent<double, Objective,
|
|
// gdc::DecreaseBacktracking<double>,
|
|
// OptimizationCallback>
|
|
// gdc::GradientDescent<double, Objective,
|
|
// gdc::DecreaseBacktracking<double>,
|
|
// OptimizationCallback>
|
|
optimizer;
|
|
|
|
// Turn verbosity on, so the optimizer prints status updates after each
|
|
// iteration.
|
|
optimizer.setVerbosity(1);
|
|
|
|
// Set initial guess.
|
|
matplot::xlabel("Optimization iterations");
|
|
matplot::ylabel("Objective value");
|
|
// matplot::figure(false);
|
|
matplot::grid(matplot::on);
|
|
// Start the optimization
|
|
auto result = optimizer.minimize(initialGuess);
|
|
|
|
std::cout << "Done! Converged: " << (result.converged ? "true" : "false")
|
|
<< " Iterations: " << result.iterations << std::endl;
|
|
|
|
// do something with final function value
|
|
std::cout << "Final fval: " << result.fval << std::endl;
|
|
|
|
// do something with final x-value
|
|
std::cout << "Final xval: " << result.xval.transpose() << std::endl;
|
|
|
|
SimulationResults reducedModelOptimizedResults =
|
|
simulator.executeSimulation(reducedModelSimulationJob);
|
|
reducedModelOptimizedResults.simulationLabel = "reducedModel";
|
|
reducedModelOptimizedResults.draw(reducedModelSimulationJob);
|
|
|
|
return result.xval;
|
|
} else { // use bobyqa
|
|
double (*pObjectiveFunction)(long, const double *) = &objective;
|
|
const size_t n = 4;
|
|
const size_t npt = 8;
|
|
assert(npt <= (n + 1) * (n + 2) / 2 && npt >= n + 2);
|
|
assert(npt < 2 * n + 1 && "The choice of the number of interpolation "
|
|
"conditions is not recommended.");
|
|
std::vector<double> x{initialGuess(0), initialGuess(1), initialGuess(2),
|
|
initialGuess(3)};
|
|
std::vector<double> xLow(x.size(), -100);
|
|
std::vector<double> xUpper(x.size(), 100);
|
|
const double maxX = *std::max_element(
|
|
x.begin(), x.end(),
|
|
[](const double &a, const double &b) { return abs(a) < abs(b); });
|
|
const double rhobeg = std::min(0.95, 0.2 * maxX);
|
|
const double rhoend = rhobeg * 1e-6;
|
|
const size_t wSize = (npt + 5) * (npt + n) + 3 * n * (n + 5) / 2;
|
|
std::vector<double> w(wSize);
|
|
bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(),
|
|
rhobeg, rhoend, 100, w.data());
|
|
|
|
std::cout << "Final objective value:" << objective(n, x.data())
|
|
<< std::endl;
|
|
|
|
Eigen::VectorXd eigenX(x.size(), 1);
|
|
|
|
for (size_t xi = 0; xi < x.size(); xi++) {
|
|
eigenX(xi) = x[xi];
|
|
}
|
|
|
|
SimulationResults reducedModelOptimizedResults =
|
|
simulator.executeSimulation(reducedModelSimulationJob);
|
|
reducedModelOptimizedResults.simulationLabel = "reducedModel";
|
|
reducedModelOptimizedResults.draw(reducedModelSimulationJob);
|
|
|
|
return eigenX;
|
|
}
|
|
}
|
|
|
|
std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
|
|
const std::shared_ptr<ElementalMesh> &pMesh) {
|
|
std::vector<SimulationJob> scenarios;
|
|
// Pull up
|
|
////1Vertex
|
|
std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
|
|
pullUp1_FixedVertices;
|
|
pullUp1_FixedVertices[3] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
|
|
std::unordered_map<VertexIndex, Vector6d> pullUp1_NodalForces{
|
|
{15, {0, 0, 100, 0, 0, 0}}};
|
|
scenarios.push_back({pMesh, pullUp1_FixedVertices, pullUp1_NodalForces, {}});
|
|
|
|
// /// 2Vertices
|
|
// std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
|
|
// pullUp1_FixedVertices;
|
|
// saddleFixedVertices[3] = std::unordered_set<DoFType>{0, 1, 2};
|
|
// saddleFixedVertices[7] = std::unordered_set<DoFType>{0, 1, 2};
|
|
// std::unordered_map<VertexIndex, Vector6d> pullUp1_NodalForces{
|
|
// {15, {0, 0, 0, 0, -900, 0}}, {3, {0, 0, 0, 0, 900, 0}},
|
|
// {7, {0, 0, 0, 700, 0, 0}}, {11, {0, 0, 0, 700, 0, 0}},
|
|
// {19, {0, 0, 0, -700, 0, 0}}, {23, {0, 0, 0, -700, 0, 0}}};
|
|
// scenarios.push_back({pMesh, pullUp1_FixedVertices, pullUp1_NodalForces,
|
|
// {}});
|
|
|
|
/// 3Vertices
|
|
// std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
|
|
// pullUp1_FixedVertices;
|
|
// saddleFixedVertices[3] = std::unordered_set<DoFType>{0, 1, 2};
|
|
// saddleFixedVertices[7] = std::unordered_set<DoFType>{0, 1, 2};
|
|
// saddleFixedVertices[11] = std::unordered_set<DoFType>{0, 1, 2};
|
|
// saddleFixedVertices[15] = std::unordered_set<DoFType>{0, 1, 2};
|
|
// saddleFixedVertices[19] = std::unordered_set<DoFType>{0, 1, 2};
|
|
// saddleFixedVertices[23] = std::unordered_set<DoFType>{0, 1, 2};
|
|
// std::unordered_map<VertexIndex, Vector6d> pullUp1_NodalForces{
|
|
// {15, {0, 0, 0, 0, -900, 0}}, {3, {0, 0, 0, 0, 900, 0}},
|
|
// {7, {0, 0, 0, 700, 0, 0}}, {11, {0, 0, 0, 700, 0, 0}},
|
|
// {19, {0, 0, 0, -700, 0, 0}}, {23, {0, 0, 0, -700, 0, 0}}};
|
|
// scenarios.push_back({pMesh, pullUp1_FixedVertices, pullUp1_NodalForces,
|
|
// {}});
|
|
|
|
// Single moment
|
|
// std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
|
|
// fixedVertices; fixedVertices[3] = std::unordered_set<DoFType>{0, 1, 2, 3,
|
|
// 4, 5}; fixedVertices[15] = std::unordered_set<DoFType>{0, 1, 2};
|
|
// std::unordered_map<VertexIndex, Vector6d> nodalForces{
|
|
// {15, {0, 0, 0, 0, 700, 0}}};
|
|
// scenarios.push_back({pMesh, saddleFixedVertices, saddleNodalForces, {}});
|
|
// scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
|
|
|
|
// Saddle
|
|
// std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
|
|
// saddle_fixedVertices;
|
|
// // saddle_fixedVertices[3] = std::unordered_set<DoFType>{0, 1, 2};
|
|
// saddle_fixedVertices[7] = std::unordered_set<DoFType>{0, 1, 2};
|
|
// saddle_fixedVertices[11] = std::unordered_set<DoFType>{0, 1, 2};
|
|
// // saddle_fixedVertices[15] = std::unordered_set<DoFType>{0, 1, 2};
|
|
// // saddle_fixedVertices[19] = std::unordered_set<DoFType>{0, 1, 2};
|
|
// // saddle_fixedVertices[23] = std::unordered_set<DoFType>{0, 1, 2};
|
|
// std::unordered_map<VertexIndex, Vector6d> saddle_nodalForces{
|
|
// {15, {0, 0, 0, 0, -1400, 0}}, {3, {0, 0, 0, 0, 1400, 0}},
|
|
// {7, {0, 0, 0, 700, 0, 0}}, {11, {0, 0, 0, 700, 0, 0}},
|
|
// {19, {0, 0, 0, -700, 0, 0}}, {23, {0, 0, 0, -700, 0, 0}}};
|
|
// // scenarios.push_back({pMesh, saddleFixedVertices, saddleNodalForces,
|
|
// {}}); std::unordered_map<VertexIndex, Eigen::Vector3d>
|
|
// saddle_forcedDisplacements; scenarios.push_back({pMesh,
|
|
// saddle_fixedVertices, saddle_nodalForces,
|
|
// saddle_forcedDisplacements});
|
|
|
|
return scenarios;
|
|
}
|
|
|
|
Eigen::VectorXd ReducedModelOptimizer::optimize() {
|
|
std::vector<SimulationJob> simulationJobs =
|
|
createScenarios(pFullModelElementalMesh);
|
|
std::vector<Eigen::VectorXd> results;
|
|
for (const SimulationJob &job : simulationJobs) {
|
|
// fullModelSimulationJob.draw();
|
|
auto result = optimizeForSimulationJob(job);
|
|
results.push_back(result);
|
|
}
|
|
|
|
if (results.empty()) {
|
|
return Eigen::VectorXd();
|
|
}
|
|
|
|
return results[0];
|
|
}
|