Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
51 changes: 24 additions & 27 deletions examples/AdjointSensitivity/AdjointSensitivity.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -57,14 +57,13 @@
*
*/


#include <iostream>
#include <iomanip>
#include <iostream>

#include <Model/PowerFlow/Bus/BusSlack.hpp>
#include <Model/PowerFlow/Generator4/Generator4.hpp>
#include <SystemModel.hpp>
#include <Solver/Dynamic/Ida.hpp>
#include <SystemModel.hpp>
#include <Utilities/Testing.hpp>

/*
Expand Down Expand Up @@ -112,7 +111,6 @@ int main()
idas->configureQuadrature();
idas->initializeQuadrature();


idas->runSimulation(0.1, 2);
idas->saveInitialCondition();

Expand All @@ -138,30 +136,29 @@ int main()
std::cout << "\n\nCost of computing objective function:\n\n";
idas->printFinalStats();

const double g1 = Q[0];
const double g1 = Q[0];
const double eps = 2e-3;

// Compute gradient of the objective function numerically
std::vector<double> dGdp(model->size_opt());

for (unsigned i=0; i<model->size_opt(); ++i)
for (unsigned i = 0; i < model->size_opt(); ++i)
{
model->param()[i] += eps;
idas->getSavedInitialCondition();
idas->initializeSimulation(t_init);
idas->initializeQuadrature();
idas->runSimulationQuadrature(t_final,100);

std::cout << "\n\nCost of computing derivative with respect to parameter "
<< i << ":\n\n";
idas->printFinalStats();
double g2 = Q[0];

// restore parameter to original value
model->param()[i] -= eps;

// Evaluate dG/dp numerically
dGdp[i] = (g2 - g1)/eps;
model->param()[i] += eps;
idas->getSavedInitialCondition();
idas->initializeSimulation(t_init);
idas->initializeQuadrature();
idas->runSimulationQuadrature(t_final, 100);

std::cout << "\n\nCost of computing derivative with respect to parameter " << i << ":\n\n";
idas->printFinalStats();
double g2 = Q[0];

// restore parameter to original value
model->param()[i] -= eps;

// Evaluate dG/dp numerically
dGdp[i] = (g2 - g1) / eps;
}

// Compute gradient of the objective function using adjoint method
Expand All @@ -185,15 +182,15 @@ int main()
int retval = 0;
std::cout << "\n\nComparison of numerical and adjoint results:\n\n";
double* neg_dGdp = idas->getAdjointIntegral();
for (unsigned i=0; i<model->size_opt(); ++i)
for (unsigned i = 0; i < model->size_opt(); ++i)
{
std::cout << "dG/dp" << i << " (numerical) = " << dGdp[i] << "\n";
std::cout << "dG/dp" << i << " (numerical) = " << dGdp[i] << "\n";
std::cout << "dG/dp" << i << " (adjoint) = " << -neg_dGdp[i] << "\n\n";
if(!isEqual(dGdp[i], -neg_dGdp[i], 10*eps))
--retval;
if (!isEqual(dGdp[i], -neg_dGdp[i], 10 * eps))
--retval;
}

if(retval < 0)
if (retval < 0)
{
std::cout << "The two results differ beyond solver tolerance!\n";
}
Expand Down
126 changes: 63 additions & 63 deletions examples/DistributedGeneratorTest/DGTest.cpp
Original file line number Diff line number Diff line change
@@ -1,85 +1,85 @@


#include <iostream>
#include <iomanip>
#include <cmath>
#include <fstream>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <math.h>

#include <Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp>


/**
* @brief Testing for the Distributed Generators outputs
*
* @param argc
* @param argv
* @return int
*
* @param argc
* @param argv
* @return int
*/
int main(int argc, char const *argv[])
int main(int argc, char const* argv[])
{

ModelLib::DistributedGeneratorParameters<double, size_t> parms;
//Parameters from MATLAB Microgrid code for first DG
parms.wb_ = 2.0*M_PI*50.0;
parms.wc_ = 31.41;
parms.mp_ = 9.4e-5;
parms.Vn_ = 380;
parms.nq_ = 1.3e-3;
parms.F_ = 0.75;
parms.Kiv_ = 420.0;
parms.Kpv_ = 0.1;
parms.Kic_ = 20.0 * 1.0e3;
parms.Kpc_ = 15.0;
parms.Cf_ = 50.0e-6;
parms.rLf_ = 0.1;
parms.Lf_ = 1.35e-3;
parms.rLc_ = 0.03;
parms.Lc_ = 0.35e-3;

ModelLib::DistributedGenerator<double, size_t> *dg = new ModelLib::DistributedGenerator<double, size_t>(0, parms, true);
ModelLib::DistributedGeneratorParameters<double, size_t> parms;
// Parameters from MATLAB Microgrid code for first DG
parms.wb_ = 2.0 * M_PI * 50.0;
parms.wc_ = 31.41;
parms.mp_ = 9.4e-5;
parms.Vn_ = 380;
parms.nq_ = 1.3e-3;
parms.F_ = 0.75;
parms.Kiv_ = 420.0;
parms.Kpv_ = 0.1;
parms.Kic_ = 20.0 * 1.0e3;
parms.Kpc_ = 15.0;
parms.Cf_ = 50.0e-6;
parms.rLf_ = 0.1;
parms.Lf_ = 1.35e-3;
parms.rLc_ = 0.03;
parms.Lc_ = 0.35e-3;

ModelLib::DistributedGenerator<double, size_t>* dg =
new ModelLib::DistributedGenerator<double, size_t>(0, parms, true);

std::vector<double> t1(16,0.0);
std::vector<double> t2{0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5};
std::vector<double> t1(16, 0.0);
std::vector<double> t2{0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5};

dg->allocate();
dg->allocate();

dg->y() = t2;
dg->yp() = t1;
dg->y() = t2;
dg->yp() = t1;

dg->evaluateResidual();
dg->evaluateResidual();

std::cout << "Output: {";
for (double i : dg->getResidual())
{
printf("%e ,", i);
}
std::cout << "}\n";
std::cout << "Output: {";
for (double i : dg->getResidual())
{
printf("%e ,", i);
}
std::cout << "}\n";

//Generated from matlab code with same parameters and inputs
std::vector<double> true_vec{3.141592277589793e+02,
8.941907747838389e-01,
1.846733023014284e+00,
3.141592277589793e+02,
1.014543000000000e+02,
-1.507680000000000e+01,
3.787993500000000e+02,
-1.300000000000000e+00,
2.899095146477517e+02,
2.939138495559215e+02,
1.507210571826699e+07,
1.659799832843673e+07,
-7.591593003913325e+03,
-8.376991073310774e+03,
3.337988298081817e+03,
2.684419146397466e+03};
// Generated from matlab code with same parameters and inputs
std::vector<double> true_vec{3.141592277589793e+02,
8.941907747838389e-01,
1.846733023014284e+00,
3.141592277589793e+02,
1.014543000000000e+02,
-1.507680000000000e+01,
3.787993500000000e+02,
-1.300000000000000e+00,
2.899095146477517e+02,
2.939138495559215e+02,
1.507210571826699e+07,
1.659799832843673e+07,
-7.591593003913325e+03,
-8.376991073310774e+03,
3.337988298081817e+03,
2.684419146397466e+03};

std::cout << "Test the Relative Error\n";
for (size_t i = 0; i < true_vec.size(); i++)
{
printf("%e ,\n", (true_vec[i] - dg->getResidual()[i]) / true_vec[i]);
}
std::cout << "Test the Relative Error\n";
for (size_t i = 0; i < true_vec.size(); i++)
{
printf("%e ,\n", (true_vec[i] - dg->getResidual()[i]) / true_vec[i]);
}

return 0;
return 0;
}
57 changes: 29 additions & 28 deletions examples/DynamicConstrainedOpt/DynamicConstrainedOpt.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -57,19 +57,17 @@
*
*/


#include <iostream>
#include <iomanip>
#include <iostream>

#include <IpIpoptApplication.hpp>
#include <IpSolveStatistics.hpp>
#include <Model/PowerFlow/Bus/BusSlack.hpp>
#include <Model/PowerFlow/Generator2/Generator2.hpp>
#include <SystemModel.hpp>
#include <Solver/Dynamic/Ida.hpp>

#include <IpIpoptApplication.hpp>
#include <IpSolveStatistics.hpp>
#include <Solver/Optimization/DynamicObjective.hpp>
#include <Solver/Optimization/DynamicConstraint.hpp>
#include <Solver/Optimization/DynamicObjective.hpp>
#include <SystemModel.hpp>
#include <Utilities/Testing.hpp>

int main()
Expand Down Expand Up @@ -114,7 +112,7 @@ int main()
{
gen->V() = 0.0;
idas->runSimulation(t_clear, 2);
gen->V() = 1.0;
gen->V() = 1.0;
gen->theta() = -0.01;
idas->saveInitialCondition();
}
Expand All @@ -131,13 +129,14 @@ int main()
// Initialize the IpoptApplication and process the options
Ipopt::ApplicationReturnStatus status;
status = ipoptApp->Initialize();
if (status != Ipopt::Solve_Succeeded) {
if (status != Ipopt::Solve_Succeeded)
{
std::cout << "\n\n*** Initialization failed! ***\n\n";
return (int) status;
}

// Set solver tolerance
const double tol = 1e-4;
const double tol = 1e-4;

// Configure Ipopt application
ipoptApp->Options()->SetStringValue("hessian_approximation", "limited-memory");
Expand All @@ -155,18 +154,19 @@ int main()
status = ipoptApp->OptimizeTNLP(ipoptDynamicObjectiveInterface);
std::cout << "\n\nProblem formulated as dynamic objective optimization ...\n";

if (status == Ipopt::Solve_Succeeded) {
if (status == Ipopt::Solve_Succeeded)
{
// Print result
std::cout << "\nSucess:\n The problem solved in "
<< ipoptApp->Statistics()->IterationCount() << " iterations!\n"
std::cout << "\nSucess:\n The problem solved in " << ipoptApp->Statistics()->IterationCount()
<< " iterations!\n"
<< " Optimal value of Pm = " << model->param()[0] << "\n"
<< " The final value of the objective function G(Pm) = "
<< ipoptApp->Statistics()->FinalObjective() << "\n\n";
<< " The final value of the objective function G(Pm) = " << ipoptApp->Statistics()->FinalObjective()
<< "\n\n";
}

// Store dynamic objective optimization results
double* results = new double[model->size_opt()];
for(unsigned i=0; i <model->size_opt(); ++i)
double* results = new double[model->size_opt()];
for (unsigned i = 0; i < model->size_opt(); ++i)
{
results[i] = model->param()[i];
}
Expand All @@ -177,34 +177,35 @@ int main()

// Initialize problem
model->param()[0] = Pm;

// Solve the problem
status = ipoptApp->OptimizeTNLP(ipoptDynamicConstraintInterface);
std::cout << "\n\nProblem formulated as dynamic constraint optimization ...\n";

if (status == Ipopt::Solve_Succeeded) {
if (status == Ipopt::Solve_Succeeded)
{
// Print result
std::cout << "\nSucess:\n The problem solved in "
<< ipoptApp->Statistics()->IterationCount() << " iterations!\n"
std::cout << "\nSucess:\n The problem solved in " << ipoptApp->Statistics()->IterationCount()
<< " iterations!\n"
<< " Optimal value of Pm = " << model->param()[0] << "\n"
<< " The final value of the objective function G(Pm) = "
<< ipoptApp->Statistics()->FinalObjective() << "\n\n";
<< " The final value of the objective function G(Pm) = " << ipoptApp->Statistics()->FinalObjective()
<< "\n\n";
}

// Compare results of the two optimization methods
int retval = 0;
for(unsigned i=0; i <model->size_opt(); ++i)
for (unsigned i = 0; i < model->size_opt(); ++i)
{
if(!isEqual(results[i], model->param()[i], 10*tol))
--retval;
if (!isEqual(results[i], model->param()[i], 10 * tol))
--retval;
}

if(retval < 0)
if (retval < 0)
{
std::cout << "The two results differ beyond solver tolerance!\n";
}

delete [] results;
delete[] results;
delete idas;
delete model;
return retval;
Expand Down
5 changes: 3 additions & 2 deletions examples/Enzyme/Library/library.cpp
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
double square(double x) {
return x * x;
double square(double x)
{
return x * x;
}
Loading
Loading