1 #ifndef ROL_DIODECIRCUIT_HPP 2 #define ROL_DIODECIRCUIT_HPP 41 typedef typename vector::size_type
uint;
47 Teuchos::RCP<std::vector<Real> >
Imeas_;
49 Teuchos::RCP<std::vector<Real> >
Vsrc_;
78 Real true_Is, Real true_Rs,
80 bool use_adjoint,
int use_hessvec)
81 : Vth_(Vth), lambertw_(lambertw), use_adjoint_(use_adjoint), use_hessvec_(use_hessvec) {
82 int n = (Vsrc_max-Vsrc_min)/Vsrc_step + 1;
83 Vsrc_ = Teuchos::rcp(
new std::vector<Real>(n,0.0));
84 Imeas_ = Teuchos::rcp(
new std::vector<Real>(n,0.0));
85 std::ofstream output (
"Measurements.dat");
86 Real left = 0.0, right = 1.0;
89 std::cout <<
"Generating data using Lambert-W function." << std::endl;
92 std::cout <<
"Generating data using Newton's method." << std::endl;
94 for (
int i = 0; i < n; i++ ) {
95 (*Vsrc_)[i] = Vsrc_min+i*Vsrc_step;
101 (*Imeas_)[i] =
Newton(I0,Vsrc_min+i*Vsrc_step,true_Is,true_Rs);
104 (*Imeas_)[i] += noise*pow(10,(
int)log10((*Imeas_)[i]))*
random(left, right);
107 if( output.is_open() ) {
108 output << std::setprecision(8) << std::scientific << (*Vsrc_)[i] <<
" " << (*Imeas_)[i] <<
"\n";
129 bool use_adjoint,
int use_hessvec)
130 : Vth_(Vth), lambertw_(lambertw), use_adjoint_(use_adjoint), use_hessvec_(use_hessvec) {
133 for(
int k = 0; std::getline(input_file,line); ++k ) {
137 input_file.seekg(0,std::ios::beg);
138 Vsrc_ = Teuchos::rcp(
new std::vector<Real>(dim,0.0));
139 Imeas_ = Teuchos::rcp(
new std::vector<Real>(dim,0.0));
141 std::cout <<
"Using input file to generate data." <<
"\n";
142 for(
int i = 0; i < dim; i++ ){
146 (*Imeas_)[i] = Imeas;
167 for ( uint i = 0; i < n; i++ ) {
175 for ( uint i = 0; i < n; i++ ) {
176 (*Ip)[i] =
Newton(I0,(*Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
189 using Teuchos::RCP;
using Teuchos::rcp;
191 uint n = Imeas_->size();
192 STDV I( rcp(
new vector(n,0.0) ) );
199 for ( uint i = 0; i < n; i++ ) {
200 val += ((*Ip)[i]-(*Imeas_)[i])*((*Ip)[i]-(*Imeas_)[i]);
208 using Teuchos::RCP;
using Teuchos::rcp;
212 uint n = Imeas_->size();
214 STDV I( rcp(
new vector(n,0.0) ) );
220 if ( use_adjoint_ ) {
222 STDV lambda( rcp(
new vector(n,0.0) ) );
229 (*gp)[0] = 0.0; (*gp)[1] = 0.0;
230 for ( uint i = 0; i < n; i++ ) {
231 (*gp)[0] +=
diodeIs((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1])*(*lambdap)[i];
232 (*gp)[1] +=
diodeRs((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1])*(*lambdap)[i];
237 STDV sensIs( rcp(
new vector(n,0.0) ) );
238 STDV sensRs( rcp(
new vector(n,0.0) ) );
247 std::ofstream output (
"Sensitivities.dat");
248 for ( uint k = 0; k < n; k++ ) {
249 if ( output.is_open() ) {
250 output << std::scientific << (*sensIsp)[k] <<
" " << (*sensRsp)[k] <<
"\n";
255 (*gp)[0] = 0.0; (*gp)[1] = 0.0;
256 for( uint i = 0; i < n; i++ ) {
257 (*gp)[0] += ((*Ip)[i]-(*Imeas_)[i])*(*sensIsp)[i];
258 (*gp)[1] += ((*Ip)[i]-(*Imeas_)[i])*(*sensRsp)[i];
272 using Teuchos::RCP;
using Teuchos::rcp;
274 if ( use_hessvec_ == 0 ) {
277 else if ( use_hessvec_ == 1 ) {
282 uint n = Imeas_->size();
284 STDV I( rcp(
new vector(n,0.0) ) );
290 STDV lambda( rcp(
new vector(n,0.0) ) );
296 STDV w( rcp(
new vector(n,0.0) ) );
300 for ( uint i = 0; i < n; i++ ){
301 (*wp)[i] = ( (*vp)[0] *
diodeIs( (*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1] )
302 + (*vp)[1] *
diodeRs( (*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1] ) )
303 /
diodeI((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
306 STDV p( rcp(
new vector(n,0.0) ) );
310 for ( uint j = 0; j < n; j++ ) {
311 (*pp)[j] = ( (*wp)[j] + (*lambdap)[j] *
diodeII( (*Ip)[j],(*Vsrc_)[j],(*Sp)[0],(*Sp)[1] )
312 * (*wp)[j] - (*lambdap)[j] *
diodeIIs( (*Ip)[j],(*Vsrc_)[j],(*Sp)[0],(*Sp)[1] )
313 * (*vp)[0] - (*lambdap)[j] *
diodeIRs( (*Ip)[j],(*Vsrc_)[j],(*Sp)[0],(*Sp)[1] )
314 * (*vp)[1] ) /
diodeI( (*Ip)[j],(*Vsrc_)[j],(*Sp)[0],(*Sp)[1] );
318 (*hvp)[0] = 0.0;(*hvp)[1] = 0.0;
319 for ( uint k = 0; k < n; k++ ) {
320 (*hvp)[0] +=
diodeIs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] )* (*pp)[k]
321 - (*lambdap)[k] * (*wp)[k] *
diodeIIs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] )
322 + (*lambdap)[k] * (*vp)[0] *
diodeIsIs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] )
323 + (*lambdap)[k] * (*vp)[1] *
diodeIsRs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] );
324 (*hvp)[1] +=
diodeRs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] ) * (*pp)[k]
325 - (*lambdap)[k] * (*wp)[k] *
diodeIRs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] )
326 + (*lambdap)[k] * (*vp)[0] *
diodeIsRs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] )
327 + (*lambdap)[k] * (*vp)[1] *
diodeRsRs( (*Ip)[k],(*Vsrc_)[k],(*Sp)[0],(*Sp)[1] );
330 else if ( use_hessvec_ == 2 ) {
336 uint n = Imeas_->size();
338 STDV I( rcp(
new vector(n,0.0) ) );
345 STDV sensIs( rcp(
new vector(n,0.0) ) );
346 STDV sensRs( rcp(
new vector(n,0.0) ) );
355 Real H11 = 0.0; Real H12 = 0.0; Real H22 = 0.0;
356 for ( uint k = 0; k < n; k++ ) {
357 H11 += (*sensIsp)[k]*(*sensIsp)[k];
358 H12 += (*sensIsp)[k]*(*sensRsp)[k];
359 H22 += (*sensRsp)[k]*(*sensRsp)[k];
363 (*hvp)[0] = H11*(*vp)[0] + H12*(*vp)[1];
364 (*hvp)[1] = H12*(*vp)[0] + H22*(*vp)[1];
382 void generate_plot(Real Is_lo, Real Is_up, Real Is_step, Real Rs_lo, Real Rs_up, Real Rs_step){
383 Teuchos::RCP<std::vector<double> > S_rcp = Teuchos::rcp(
new std::vector<double>(2,0.0) );
385 std::ofstream output (
"Objective.dat");
391 int n = (Is_up-Is_lo)/Is_step + 1;
392 int m = (Rs_up-Rs_lo)/Rs_step + 1;
393 for (
int i = 0; i < n; i++ ) {
394 Is = Is_lo + i*Is_step;
395 for (
int j = 0; j < m; j++ ) {
396 Rs = Rs_lo + j*Rs_step;
400 if ( output.is_open() ) {
401 output << std::scientific << Is <<
" " << Rs <<
" " << val << std::endl;
411 using Teuchos::dyn_cast;
using Teuchos::getConst;
413 return dyn_cast<
const STDV>(getConst(x)).
getVector();
415 catch (std::exception &e) {
417 return dyn_cast<
const PSV>(getConst(x)).
getVector();
419 catch (std::exception &e) {
420 return dyn_cast<
const DSV>(getConst(x)).
getVector();
426 using Teuchos::dyn_cast;
430 catch (std::exception &e) {
434 catch (std::exception &e) {
440 Real
random(
const Real left,
const Real right)
const {
441 return (Real)rand()/(Real)RAND_MAX * (right - left) + left;
454 Real
diode(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
455 return I-Is*(exp((Vsrc-I*Rs)/Vth_)-1);
459 Real
diodeI(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
460 return 1+Is*exp((Vsrc-I*Rs)/Vth_)*(Rs/
Vth_);
464 Real
diodeIs(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
465 return 1-exp((Vsrc-I*Rs)/Vth_);
469 Real
diodeRs(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
470 return Is*exp((Vsrc-I*Rs)/Vth_)*(I/
Vth_);
474 Real
diodeII(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
475 return -Is*exp((Vsrc-I*Rs)/Vth_)*(Rs/
Vth_)*(Rs/Vth_);
479 Real
diodeIIs(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
480 return exp((Vsrc-I*Rs)/Vth_)*(Rs/
Vth_);
484 Real
diodeIRs(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
485 return (Is/Vth_)*exp((Vsrc-I*Rs)/Vth_)*(1-(I*Rs)/Vth_);
489 Real
diodeIsIs(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
494 Real
diodeIsRs(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
495 return exp((Vsrc-I*Rs)/Vth_)*(I/
Vth_);
499 Real
diodeRsRs(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
500 return -Is*exp((Vsrc-I*Rs)/Vth_)*(I/
Vth_)*(I/Vth_);
510 Real
Newton(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
515 Real fval =
diode(IN,Vsrc,Is,Rs);
520 for (
int i = 0; i < MAXIT; i++ ) {
521 if ( std::abs(fval) < TOL ) {
525 dfval =
diodeI(IN,Vsrc,Is,Rs);
526 if( std::abs(dfval) < EPS ){
527 std::cout <<
"denominator is too small" << std::endl;
532 IN_tmp = IN - alpha*fval/dfval;
533 fval_tmp =
diode(IN_tmp,Vsrc,Is,Rs);
534 while ( std::abs(fval_tmp) >= (1.0-1.e-4*alpha)*std::abs(fval) ) {
536 IN_tmp = IN - alpha*fval/dfval;
537 fval_tmp =
diode(IN_tmp,Vsrc,Is,Rs);
538 if ( alpha < std::sqrt(EPS) ) {
583 void lambertw(Real x, Real &w,
int &ierr, Real &xi){
584 int i = 0, maxit = 10;
585 const Real turnpt = -exp(-1.), c1 = 1.5, c2 = .75;
586 Real r, r2, r3, s, mach_eps, relerr = 1., diff;
601 w = x*(1.-x + c1*x*x);
608 w = x*(1.0-x + c1*x*x);
609 xi = log(1.0-x + c1*x*x) - w;
616 w = -1 + sqrt(2.0*exp(1.))*sqrt(x-turnpt);
631 while ( relerr > mach_eps && i < maxit ) {
635 s = 6.*(w+1.0)*(w+1.0);
636 w = w * ( 1.0 + r + r2/(2.0*( w+1.0)) - (2. * w -1.0)*r3/s );
637 w = ((w*x < 0.0) ? -w : w);
640 relerr = ((x > 1.0) ? xi/w : xi);
641 relerr = ((relerr < 0.0) ? -relerr : relerr);
644 ierr = ((i == maxit) ? 2 : ierr);
660 Real arg1 = (Vsrc + Is*Rs)/Vth_;
661 Real evd = exp(arg1);
662 Real lambWArg = Is*Rs*evd/
Vth_;
663 Real lambWReturn = 0.0;
664 Real lambWError = 0.0;
666 lambertw(lambWArg, lambWReturn, ierr, lambWError);
668 std::cout <<
"LambertW error: argument is not in the domain" << std::endl;
672 std::cout <<
"LambertW error: BUG!" << std::endl;
674 Real Id = -Is+Vth_*(lambWReturn)/Rs;
694 for ( uint i = 0; i < n; i++ ){
695 (*lambdap)[i] = ((*Imeas_)[i]-(*Ip)[i])
696 /
diodeI((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
715 for ( uint i = 0; i < n; i++ ) {
716 (*sensp)[i] = -
diodeIs((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1])
717 /
diodeI((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
736 for ( uint i = 0; i < n; i++ ) {
737 (*sensp)[i] = -
diodeRs((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1])
738 /
diodeI((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
Provides the interface to evaluate objective functions.
bool use_adjoint_
If true, use adjoint gradient computation, else compute gradient using sensitivities.
Real diodeIsRs(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Second derivative of diode equation wrt Is and Rs.
Real noise_
Percentage of noise to add to measurements; if 0.0 - no noise.
Real diodeIsIs(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Second derivative of diode equation wrt Is^2.
std::vector< Real > vector
Objective_DiodeCircuit(Real Vth, std::ifstream &input_file, bool lambertw, Real noise, bool use_adjoint, int use_hessvec)
A constructor using data from given file.
Real lambertWCurrent(Real Is, Real Rs, Real Vsrc)
Find currents using Lambert-W function.
Provides the std::vector implementation of the ROL::Vector interface that handles scalings in the inn...
Teuchos::RCP< vector > getVector(V &x)
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
void solve_sensitivity_Is(Vector< Real > &sens, const Vector< Real > &I, const Vector< Real > &S)
Solve the sensitivity equation wrt Is.
void solve_adjoint(Vector< Real > &lambda, const Vector< Real > &I, const Vector< Real > &S)
Solve the adjoint equation.
void gradient(Vector< Real > &g, const Vector< Real > &S, Real &tol)
Compute the gradient of the reduced objective function either using adjoint or using sensitivities...
Teuchos::RCP< std::vector< Real > > Imeas_
Vector of measured currents in DC analysis (data)
Real Vth_
Thermal voltage (constant)
void generate_plot(Real Is_lo, Real Is_up, Real Is_step, Real Rs_lo, Real Rs_up, Real Rs_step)
Generate data to plot objective function.
Defines the linear algebra or vector space interface.
Teuchos::RCP< std::vector< Real > > Vsrc_
Vector of source voltages in DC analysis (input)
The diode circuit problem.
Real value(const Vector< Real > &S, Real &tol)
Evaluate objective function.
Real diodeI(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Derivative of diode equation wrt I.
Real random(const Real left, const Real right) const
void lambertw(Real x, Real &w, int &ierr, Real &xi)
Lambert-W function for diodes.
Real Newton(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Newton's method with line search.
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &S, Real &tol)
Compute the Hessian-vector product of the reduced objective function.
void solve_circuit(Vector< Real > &I, const Vector< Real > &S)
Solve circuit given optimization parameters Is and Rs.
Provides the std::vector implementation of the ROL::Vector interface that handles scalings in the inn...
PrimalScaledStdVector< Real > PSV
bool lambertw_
If true, use Lambert-W function to solve circuit, else use Newton's method.
Real diode(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Diode equation.
Real diodeIs(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Derivative of diode equation wrt Is.
Real diodeRsRs(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Second derivative of diode equation wrt Rs^2.
Objective_DiodeCircuit(Real Vth, Real Vsrc_min, Real Vsrc_max, Real Vsrc_step, Real true_Is, Real true_Rs, bool lambertw, Real noise, bool use_adjoint, int use_hessvec)
A constructor generating data.
void solve_sensitivity_Rs(Vector< Real > &sens, const Vector< Real > &I, const Vector< Real > &S)
Solve the sensitivity equation wrt Rs.
Real diodeIIs(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Second derivative of diode equation wrt I and Is.
Real diodeIRs(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Second derivative of diode equation wrt I and Rs.
Real diodeII(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Second derivative of diode equation wrt I^2.
Teuchos::RCP< const vector > getVector(const V &x)
DualScaledStdVector< Real > DSV
Real diodeRs(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Derivative of diode equation wrt Rs.
void set_method(bool lambertw)
Change the method for solving the circuit if needed.