Skip to content

Commit

Permalink
Merge pull request #146 from usnistgov/pureparamoptim
Browse files Browse the repository at this point in the history 
Pure Parameter Optimization
  • Loading branch information
@ianhbell
ianhbell authored Jul 23, 2024
2 parents 98f7f14 + dfdf55d commit 3112bbe
Show file tree
Hide file tree
Showing 7 changed files with 559 additions and 9 deletions.
180 changes: 180 additions & 0 deletions doc/source/recipes/FitPureParameters.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,180 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"id": "0bace5e5",
"metadata": {},
"outputs": [],
"source": [
"import teqp\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import CoolProp.CoolProp as CP\n",
"import scipy.optimize\n",
"display(teqp.__version__)\n",
"\n",
"nonpolar = {\n",
" \"kind\": \"SAFT-VR-Mie\",\n",
" \"model\": {\n",
" \"coeffs\": [\n",
" {\n",
" \"name\": \"R32\",\n",
" \"BibTeXKey\": \"Bell\",\n",
" \"m\": 1.2476268271391935,\n",
" \"sigma_m\": 3.6080717234117107e-10,\n",
" \"epsilon_over_k\": 172.53065054286867,\n",
" \"lambda_r\": 14.634722358167384,\n",
" \"lambda_a\": 6\n",
" }\n",
" ]\n",
" }\n",
"}\n",
"template = { \n",
" 'kind': 'genericSAFT', \n",
" 'model': {\n",
" 'nonpolar': nonpolar\n",
" }\n",
"}\n",
"\n",
"# NOTE: '/' inside field names MUST be escaped as ~1; see https://datatracker.ietf.org/doc/html/rfc6901#section-3\n",
"pointers = [\n",
" '/model/nonpolar/model/coeffs/0/m',\n",
" '/model/nonpolar/model/coeffs/0/sigma_m',\n",
" '/model/nonpolar/model/coeffs/0/epsilon_over_k',\n",
" '/model/nonpolar/model/coeffs/0/lambda_r',\n",
" '/model/nonpolar/model/coeffs/0/lambda_a'\n",
"]\n",
"x0 = [1.5, 3e-10, 150, 19, 5.7]\n",
"bounds = [(1,5), (2e-10,5e-10), (100,400), (12,50), (5.1, 6.0)]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f5463428",
"metadata": {},
"outputs": [],
"source": [
"FLD = 'Methane'\n",
"ppo = teqp.PureParameterOptimizer(template, pointers)\n",
"\n",
"Ts = np.linspace(100, 170, 10)\n",
"\n",
"# Generate some artificial pseudo-experimental data from \n",
"# the reference EOS as implemented in CoolProp\n",
"for T in Ts:\n",
" pt = teqp.SatRhoLPWPoint()\n",
" rhoL, rhoV = [CP.PropsSI('Dmolar','Q',Q,'T',T,FLD) for Q in [0,1]]\n",
" p = CP.PropsSI('P','Q',0,'T',T,FLD)\n",
" w = CP.PropsSI('speed_of_sound','Q',0,'T',T,FLD)\n",
" pt.T = T\n",
" # Measurands (here, pseudo-experimental values)\n",
" pt.p_exp = p\n",
" pt.rhoL_exp = rhoL\n",
" pt.w_exp = w\n",
" \n",
" # Additional parameters\n",
" pt.rhoL_guess = rhoL\n",
" pt.rhoV_guess = rhoV\n",
" pt.R = 8.31446261815324\n",
" pt.M = CP.PropsSI('molemass',FLD)\n",
" AS = CP.AbstractState('HEOS',FLD)\n",
" AS.update(CP.DmolarT_INPUTS, 1e-10, T)\n",
" pt.Ao20 = AS.tau()**2*AS.d2alpha0_dTau2() # -cv0/R\n",
" \n",
" # Weights (multiplied by 100 to put on a percentage basis)\n",
" pt.weight_p = 1*100\n",
" pt.weight_rho = 1*100\n",
" pt.weight_w = 0.25*100\n",
" ppo.add_one_contribution(pt)\n",
"\n",
"def cost_function(x):\n",
"# return ppo.cost_function_threaded(x, 10) # This is an option if you have lots of threads\n",
" return ppo.cost_function(x)\n",
"\n",
"r = scipy.optimize.differential_evolution(cost_function, bounds=bounds, disp=True, maxiter=10000, popsize=8)\n",
"print(r)\n",
"x = r.x\n",
"model = teqp.make_model(ppo.build_JSON(x))\n",
"\n",
"Tc, rhoc = model.solve_pure_critical(400, 5000)\n",
"print(Tc, rhoc)\n",
"anc = teqp.build_ancillaries(model, Tc, rhoc, 0.5*Tc)\n",
"\n",
"def _get_SOS(model, T, rho, z, *, R, M, Ao20):\n",
" \"\"\" Helper function to calculate speed of sound \"\"\"\n",
" Ar0n = model.get_Ar02n(T, rho, z)\n",
" Ar01 = Ar0n[1]; Ar02 = Ar0n[2]\n",
" Ar11 = model.get_Ar11(T, rho, z)\n",
" Ar20 = model.get_Ar20(T, rho, z)\n",
"\n",
" # M*w^2/(R*T) where w is the speed of sound\n",
" # from the definition w = sqrt(dp/drho|s)\n",
" Mw2RT = 1 + 2*Ar01 + Ar02 - (1 + Ar01 - Ar11)**2/(Ao20 + Ar20)\n",
" if Mw2RT < 0:\n",
" return 1e6\n",
" w = (Mw2RT*R*T/M)**0.5\n",
" return w\n",
"\n",
"TcREF = CP.PropsSI('Tcrit', FLD)\n",
"Tsverify = np.linspace(100, TcREF*0.99999, 1000)\n",
"RHOL, RHOV, PPP, WWW = [],[],[],[]\n",
"z = np.array([1.0])\n",
"for T in Tsverify:\n",
" rhoL, rhoV = model.pure_VLE_T(T, anc.rhoL(T)*1.01, anc.rhoV(T)*0.9, 10)\n",
" pL = rhoL*model.get_R(z)*T*(1+model.get_Ar01(T, rhoL, z))\n",
" pV = rhoV*model.get_R(z)*T*(1+model.get_Ar01(T, rhoV, z))\n",
" RHOL.append(rhoL)\n",
" RHOV.append(rhoV)\n",
" PPP.append(pL)\n",
" AS = CP.AbstractState('HEOS',FLD)\n",
" AS.update(CP.DmolarT_INPUTS, 1e-10, T)\n",
" Ao20 = AS.d2alpha0_dTau2()*AS.tau()**2\n",
" WWW.append(_get_SOS(model, T, rhoL, z, R=8.314462618,M=CP.PropsSI('molemass',FLD),Ao20=Ao20))\n",
" \n",
"# Plot the T-rho for VLE\n",
"line, = plt.plot(np.array(RHOL), Tsverify)\n",
"plt.plot(np.array(RHOV), Tsverify, color=line.get_color())\n",
"for Q in [0, 1]:\n",
" D = CP.PropsSI('Dmolar','T',Tsverify,'Q',Q,FLD)\n",
" plt.plot(D, Tsverify, lw=2, color='k')\n",
"plt.gca().set(xlabel=r'$\\rho$ / mol/m$^3$', ylabel='$T$ / K')\n",
"\n",
"# And a deviation plot, much closer to the critical point\n",
"# than the fitting region\n",
"fig, (ax1, ax2, ax3) = plt.subplots(3,1, figsize=(10,5), sharex=True)\n",
"\n",
"ax1.plot(Tsverify, (np.array(PPP)/CP.PropsSI('P','T',Tsverify,'Q',0,FLD)-1)*100)\n",
"ax1.set(ylabel=r'$(p_{fit}/p_{\\rm pexp}-1)\\times 100$')\n",
"\n",
"ax2.plot(Tsverify, (np.array(RHOL)/CP.PropsSI('Dmolar','T',Tsverify,'Q',0,FLD)-1)*100)\n",
"ax2.set(ylabel=r'$(\\rho_{fit}/\\rho_{\\rm pexp}-1)\\times 100$')\n",
"\n",
"ax3.plot(Tsverify, (np.array(WWW)/CP.PropsSI('speed_of_sound','T',Tsverify,'Q',0,FLD)-1)*100)\n",
"ax3.set(ylabel=r'$(w_{fit}/w_{\\rm pexp}-1)\\times 100$', xlabel='$T$ / K')"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
3 changes: 2 additions & 1 deletion doc/source/recipes/index.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,4 +8,5 @@ The examples here demonstrate various capabilities of teqp to solve practical pr
:caption: Contents:

HeatPumpModel

IdealCurves
FitPureParameters
196 changes: 196 additions & 0 deletions include/teqp/algorithms/pure_param_optimization.hpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,196 @@
#pragma once

#include "nlohmann/json.hpp"
#include "teqp/cpp/teqpcpp.hpp"
#include "teqp/exceptions.hpp"
#include "teqp/math/pow_templates.hpp"

#include <boost/asio/thread_pool.hpp>
#include <boost/asio/post.hpp>

#include <tuple>

namespace teqp::algorithms::pure_param_optimization {

struct SatRhoLPoint{
double T, rhoL_exp, rhoL_guess, rhoV_guess, weight=1.0;
auto check_fields() const{}
template<typename Model>
auto calculate_contribution(const Model& model) const{
auto rhoLrhoV = model->pure_VLE_T(T, rhoL_guess, rhoV_guess, 10);
return std::abs(rhoLrhoV[0]-rhoL_exp)*weight;
}
};

#define stdstringify(s) std::string(#s)

#define SatRhoLPPoint_optionalfields X(T) X(p_exp) X(rhoL_exp) X(rhoL_guess) X(rhoV_guess)
struct SatRhoLPPoint{
#define X(field) std::optional<double> field;
SatRhoLPPoint_optionalfields
#undef X
double weight_rho=1.0, weight_p=1.0, R=8.31446261815324;
Eigen::ArrayXd z = (Eigen::ArrayXd(1) << 1.0).finished();

auto check_fields() const{
#define X(field) if (!field){ throw teqp::InvalidArgument("A field [" + stdstringify(field) + "] has not been initialized"); }
SatRhoLPPoint_optionalfields
#undef X
}

template<typename Model>
auto calculate_contribution(const Model& model) const{
auto rhoLrhoV = model->pure_VLE_T(T.value(), rhoL_guess.value(), rhoV_guess.value(), 10);
auto rhoL = rhoLrhoV[0];
auto p = rhoL*R*T.value()*(1+model->get_Ar01(T.value(), rhoL, z));
// std::cout << p << "," << p_exp << "," << (p_exp-p)/p_exp << std::endl;
return std::abs(rhoL-rhoL_exp.value())/rhoL_exp.value()*weight_rho + std::abs(p-p_exp.value())/p_exp.value()*weight_p;
}
};

#define SatRhoLPWPoint_optionalfields X(T) X(p_exp) X(rhoL_exp) X(w_exp) X(rhoL_guess) X(rhoV_guess) X(Ao20) X(M) X(R)
struct SatRhoLPWPoint{

#define X(field) std::optional<double> field;
SatRhoLPWPoint_optionalfields
#undef X
double weight_rho=1.0, weight_p=1.0, weight_w=1.0;
Eigen::ArrayXd z = (Eigen::ArrayXd(1) << 1.0).finished();

auto check_fields() const{
#define X(field) if (!field){ throw teqp::InvalidArgument("A field [" + stdstringify(field) + "] has not been initialized"); }
SatRhoLPWPoint_optionalfields
#undef X
}

template<typename Model>
auto calculate_contribution(const Model& model) const{

auto rhoLrhoV = model->pure_VLE_T(T.value(), rhoL_guess.value(), rhoV_guess.value(), 10);

// First part, density
auto rhoL = rhoLrhoV[0];

auto Ar0n = model->get_Ar02n(T.value(), rhoL, z);
double Ar01 = Ar0n[1], Ar02 = Ar0n[2];
auto Ar11 = model->get_Ar11(T.value(), rhoL, z);
auto Ar20 = model->get_Ar20(T.value(), rhoL, z);

// Second part, presure
auto p = rhoL*R.value()*T.value()*(1+Ar01);

// Third part, speed of sound
//
// M*w^2/(R*T) where w is the speed of sound
// from the definition w = sqrt(dp/drho|s)
double Mw2RT = 1 + 2*Ar01 + Ar02 - POW2(1.0 + Ar01 - Ar11)/(Ao20.value() + Ar20);
double w = sqrt(Mw2RT*R.value()*T.value()/M.value());

// std::cout << p << "," << p_exp << "," << (p_exp-p)/p_exp << std::endl;
double cost_rhoL = std::abs(rhoL-rhoL_exp.value())/rhoL_exp.value()*weight_rho;
double cost_p = std::abs(p-p_exp.value())/p_exp.value()*weight_p;
double cost_w = std::abs(w-w_exp.value())/w_exp.value()*weight_w;
return cost_rhoL + cost_p + cost_w;
}
};
using PureOptimizationContribution = std::variant<SatRhoLPoint, SatRhoLPPoint, SatRhoLPWPoint>;

class PureParameterOptimizer{
private:
auto make_pointers(const std::vector<std::variant<std::string, std::vector<std::string>>>& pointerstrs){
std::vector<std::vector<nlohmann::json::json_pointer>> pointers_;
for (auto & pointer: pointerstrs){
if (std::holds_alternative<std::string>(pointer)){
const std::string& s= std::get<std::string>(pointer);
pointers_.emplace_back(1, nlohmann::json::json_pointer(s));
}
else{
std::vector<nlohmann::json::json_pointer> buffer;
for (const auto& s : std::get<std::vector<std::string>>(pointer)){
buffer.emplace_back(s);
}
pointers_.push_back(buffer);
}
}
return pointers_;
}
public:
const nlohmann::json jbase;
std::vector<std::vector<nlohmann::json::json_pointer>> pointers;
std::vector<PureOptimizationContribution> contributions;

PureParameterOptimizer(const nlohmann::json jbase, const std::vector<std::variant<std::string, std::vector<std::string>>>& pointerstrs) : jbase(jbase), pointers(make_pointers(pointerstrs)){}

void add_one_contribution(const PureOptimizationContribution& cont){
std::visit([](const auto&c){c.check_fields();}, cont);
contributions.push_back(cont);
}

auto prepare_helper_models(const auto& model) const {
return std::vector<double>{1, 3};
}

template<typename T>
nlohmann::json build_JSON(const T& x) const{

if (x.size() != pointers.size()){
throw teqp::InvalidArgument("sizes don't match");
};
nlohmann::json j = jbase;
for (auto i = 0; i < x.size(); ++i){
for (const auto& ptr : pointers[i]){
j.at(ptr) = x[i];
}
}
return j;
}

template<typename T>
auto prepare(const T& x) const {
auto j = build_JSON(x);
auto model = teqp::cppinterface::make_model(j);
auto helpers = prepare_helper_models(model);
return std::make_tuple(std::move(model), helpers);
}

template<typename T>
auto cost_function(const T& x) const{
auto [model, helpers] = prepare(x);
double cost = 0.0;
for (const auto& contrib : contributions){
cost += std::visit([&model](const auto& c){ return c.calculate_contribution(model); }, contrib);
}
if (!std::isfinite(cost)){
return 1e30;
}
return cost;
}

template<typename T>
auto cost_function_threaded(const T& x, std::size_t Nthreads) {
boost::asio::thread_pool pool{Nthreads}; // Nthreads in the pool
const auto [model, helpers] = prepare(x);
std::valarray<double> buffer(contributions.size());
std::size_t i = 0;
for (const auto& contrib : contributions){
auto& dest = buffer[i];
auto payload = [&model, &dest, contrib] (){
dest = std::visit([&model](const auto& c){ return c.calculate_contribution(model); }, contrib);
if (!std::isfinite(dest)){ dest = 1e30; }
};
boost::asio::post(pool, payload);
i++;
}
pool.join();
double summer = 0.0;
for (auto i = 0; i < contributions.size(); ++i){
// std::cout << buffer[i] << std::endl;
summer += buffer[i];
}
// std::cout << summer << std::endl;
return summer;
}

};

}
Loading

0 comments on commit 3112bbe

Please sign in to comment.