1 #ifndef ROL_DIODECIRCUIT_HPP
2 #define ROL_DIODECIRCUIT_HPP
48 ROL::Ptr<std::vector<Real> >
Vsrc_;
77 Real true_Is, Real true_Rs,
79 bool use_adjoint,
int use_hessvec)
81 int n = (Vsrc_max-Vsrc_min)/Vsrc_step + 1;
82 Vsrc_ = ROL::makePtr<std::vector<Real>>(n,0.0);
83 Imeas_ = ROL::makePtr<std::vector<Real>>(n,0.0);
84 std::ofstream output (
"Measurements.dat");
85 Real left = 0.0, right = 1.0;
88 std::cout <<
"Generating data using Lambert-W function." << std::endl;
91 std::cout <<
"Generating data using Newton's method." << std::endl;
93 for (
int i = 0; i < n; i++ ) {
94 (*Vsrc_)[i] = Vsrc_min+i*Vsrc_step;
100 (*Imeas_)[i] =
Newton(I0,Vsrc_min+i*Vsrc_step,true_Is,true_Rs);
103 (*Imeas_)[i] += noise*pow(10,(
int)log10((*
Imeas_)[i]))*
random(left, right);
106 if( output.is_open() ) {
107 output << std::setprecision(8) << std::scientific << (*Vsrc_)[i] <<
" " << (*Imeas_)[i] <<
"\n";
128 bool use_adjoint,
int use_hessvec)
132 for(
int k = 0; std::getline(input_file,line); ++k ) {
136 input_file.seekg(0,std::ios::beg);
137 Vsrc_ = ROL::makePtr<std::vector<Real>>(
dim,0.0);
138 Imeas_ = ROL::makePtr<std::vector<Real>>(
dim,0.0);
140 std::cout <<
"Using input file to generate data." <<
"\n";
141 for(
int i = 0; i <
dim; i++ ){
145 (*Imeas_)[i] = Imeas;
159 ROL::Ptr<const vector> Sp =
getVector(S);
166 for (
uint i = 0; i < n; i++ ) {
174 for (
uint i = 0; i < n; i++ ) {
175 (*Ip)[i] =
Newton(I0,(*
Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
189 ROL::Ptr<const vector> Sp =
getVector(S);
191 STDV I( ROL::makePtr<vector>(n,0.0) );
198 for (
uint i = 0; i < n; i++ ) {
199 val += ((*Ip)[i]-(*Imeas_)[i])*((*Ip)[i]-(*Imeas_)[i]);
209 ROL::Ptr<const vector> Sp =
getVector(S);
213 STDV I( ROL::makePtr<vector>(n,0.0) );
221 STDV lambda( ROL::makePtr<vector>(n,0.0) );
222 ROL::Ptr<vector> lambdap =
getVector(lambda);
228 (*gp)[0] = 0.0; (*gp)[1] = 0.0;
229 for (
uint i = 0; i < n; i++ ) {
230 (*gp)[0] +=
diodeIs((*Ip)[i],(*
Vsrc_)[i],(*Sp)[0],(*Sp)[1])*(*lambdap)[i];
231 (*gp)[1] +=
diodeRs((*Ip)[i],(*
Vsrc_)[i],(*Sp)[0],(*Sp)[1])*(*lambdap)[i];
236 STDV sensIs( ROL::makePtr<vector>(n,0.0) );
237 STDV sensRs( ROL::makePtr<vector>(n,0.0) );
242 ROL::Ptr<vector> sensIsp =
getVector(sensIs);
243 ROL::Ptr<vector> sensRsp =
getVector(sensRs);
246 std::ofstream output (
"Sensitivities.dat");
247 for (
uint k = 0; k < n; k++ ) {
248 if ( output.is_open() ) {
249 output << std::scientific << (*sensIsp)[k] <<
" " << (*sensRsp)[k] <<
"\n";
254 (*gp)[0] = 0.0; (*gp)[1] = 0.0;
255 for(
uint i = 0; i < n; i++ ) {
256 (*gp)[0] += ((*Ip)[i]-(*Imeas_)[i])*(*sensIsp)[i];
257 (*gp)[1] += ((*Ip)[i]-(*Imeas_)[i])*(*sensRsp)[i];
278 ROL::Ptr<const vector> vp =
getVector(v);
279 ROL::Ptr<const vector> Sp =
getVector(S);
283 STDV I( ROL::makePtr<vector>(n,0.0) );
289 STDV lambda( ROL::makePtr<vector>(n,0.0) );
290 ROL::Ptr<vector> lambdap =
getVector(lambda);
295 STDV w( ROL::makePtr<vector>(n,0.0) );
299 for (
uint i = 0; i < n; i++ ){
300 (*wp)[i] = ( (*vp)[0] *
diodeIs( (*Ip)[i],(*
Vsrc_)[i],(*Sp)[0],(*Sp)[1] )
301 + (*vp)[1] *
diodeRs( (*Ip)[i],(*
Vsrc_)[i],(*Sp)[0],(*Sp)[1] ) )
302 /
diodeI((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
305 STDV p( ROL::makePtr<vector>(n,0.0) );
309 for (
uint j = 0; j < n; j++ ) {
310 (*pp)[j] = ( (*wp)[j] + (*lambdap)[j] *
diodeII( (*Ip)[j],(*
Vsrc_)[j],(*Sp)[0],(*Sp)[1] )
311 * (*wp)[j] - (*lambdap)[j] *
diodeIIs( (*Ip)[j],(*
Vsrc_)[j],(*Sp)[0],(*Sp)[1] )
312 * (*vp)[0] - (*lambdap)[j] *
diodeIRs( (*Ip)[j],(*
Vsrc_)[j],(*Sp)[0],(*Sp)[1] )
313 * (*vp)[1] ) /
diodeI( (*Ip)[j],(*Vsrc_)[j],(*Sp)[0],(*Sp)[1] );
317 (*hvp)[0] = 0.0;(*hvp)[1] = 0.0;
318 for (
uint k = 0; k < n; k++ ) {
319 (*hvp)[0] +=
diodeIs( (*Ip)[k],(*
Vsrc_)[k],(*Sp)[0],(*Sp)[1] )* (*pp)[k]
320 - (*lambdap)[k] * (*wp)[k] *
diodeIIs( (*Ip)[k],(*
Vsrc_)[k],(*Sp)[0],(*Sp)[1] )
321 + (*lambdap)[k] * (*vp)[0] *
diodeIsIs( (*Ip)[k],(*
Vsrc_)[k],(*Sp)[0],(*Sp)[1] )
322 + (*lambdap)[k] * (*vp)[1] *
diodeIsRs( (*Ip)[k],(*
Vsrc_)[k],(*Sp)[0],(*Sp)[1] );
323 (*hvp)[1] +=
diodeRs( (*Ip)[k],(*
Vsrc_)[k],(*Sp)[0],(*Sp)[1] ) * (*pp)[k]
324 - (*lambdap)[k] * (*wp)[k] *
diodeIRs( (*Ip)[k],(*
Vsrc_)[k],(*Sp)[0],(*Sp)[1] )
325 + (*lambdap)[k] * (*vp)[0] *
diodeIsRs( (*Ip)[k],(*
Vsrc_)[k],(*Sp)[0],(*Sp)[1] )
326 + (*lambdap)[k] * (*vp)[1] *
diodeRsRs( (*Ip)[k],(*
Vsrc_)[k],(*Sp)[0],(*Sp)[1] );
332 ROL::Ptr<const vector> vp =
getVector(v);
333 ROL::Ptr<const vector> Sp =
getVector(S);
337 STDV I( ROL::makePtr<vector>(n,0.0) );
344 STDV sensIs( ROL::makePtr<vector>(n,0.0) );
345 STDV sensRs( ROL::makePtr<vector>(n,0.0) );
350 ROL::Ptr<vector> sensIsp =
getVector(sensIs);
351 ROL::Ptr<vector> sensRsp =
getVector(sensRs);
354 Real H11 = 0.0; Real H12 = 0.0; Real H22 = 0.0;
355 for (
uint k = 0; k < n; k++ ) {
356 H11 += (*sensIsp)[k]*(*sensIsp)[k];
357 H12 += (*sensIsp)[k]*(*sensRsp)[k];
358 H22 += (*sensRsp)[k]*(*sensRsp)[k];
362 (*hvp)[0] = H11*(*vp)[0] + H12*(*vp)[1];
363 (*hvp)[1] = H12*(*vp)[0] + H22*(*vp)[1];
381 void generate_plot(Real Is_lo, Real Is_up, Real Is_step, Real Rs_lo, Real Rs_up, Real Rs_step){
382 ROL::Ptr<std::vector<Real> > S_ptr = ROL::makePtr<std::vector<Real>>(2,0.0);
384 std::ofstream output (
"Objective.dat");
390 int n = (Is_up-Is_lo)/Is_step + 1;
391 int m = (Rs_up-Rs_lo)/Rs_step + 1;
392 for (
int i = 0; i < n; i++ ) {
393 Is = Is_lo + i*Is_step;
394 for (
int j = 0; j < m; j++ ) {
395 Rs = Rs_lo + j*Rs_step;
399 if ( output.is_open() ) {
400 output << std::scientific << Is <<
" " << Rs <<
" " << val << std::endl;
413 catch (std::exception &e) {
417 catch (std::exception &e) {
428 catch (std::exception &e) {
432 catch (std::exception &e) {
438 Real
random(
const Real left,
const Real right)
const {
439 return (Real)rand()/(Real)RAND_MAX * (right - left) + left;
452 Real
diode(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
453 return I-Is*(exp((Vsrc-I*Rs)/
Vth_)-1);
457 Real
diodeI(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
458 return 1+Is*exp((Vsrc-I*Rs)/
Vth_)*(Rs/
Vth_);
462 Real
diodeIs(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
463 return 1-exp((Vsrc-I*Rs)/
Vth_);
467 Real
diodeRs(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
468 return Is*exp((Vsrc-I*Rs)/
Vth_)*(I/
Vth_);
472 Real
diodeII(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
477 Real
diodeIIs(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
478 return exp((Vsrc-I*Rs)/
Vth_)*(Rs/
Vth_);
482 Real
diodeIRs(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
487 Real
diodeIsIs(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
492 Real
diodeIsRs(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
493 return exp((Vsrc-I*Rs)/
Vth_)*(I/
Vth_);
497 Real
diodeRsRs(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
508 Real
Newton(
const Real I,
const Real Vsrc,
const Real Is,
const Real Rs){
513 Real fval =
diode(IN,Vsrc,Is,Rs);
518 for (
int i = 0; i < MAXIT; i++ ) {
519 if ( std::abs(fval) < TOL ) {
523 dfval =
diodeI(IN,Vsrc,Is,Rs);
524 if( std::abs(dfval) < EPS ){
525 std::cout <<
"denominator is too small" << std::endl;
530 IN_tmp = IN - alpha*fval/dfval;
531 fval_tmp =
diode(IN_tmp,Vsrc,Is,Rs);
532 while ( std::abs(fval_tmp) >= (1.0-1.e-4*alpha)*std::abs(fval) ) {
534 IN_tmp = IN - alpha*fval/dfval;
535 fval_tmp =
diode(IN_tmp,Vsrc,Is,Rs);
536 if ( alpha < std::sqrt(EPS) ) {
581 void lambertw(Real x, Real &w,
int &ierr, Real &xi){
582 int i = 0, maxit = 10;
583 const Real turnpt = -exp(-1.), c1 = 1.5, c2 = .75;
584 Real r, r2, r3, s, mach_eps, relerr = 1., diff;
599 w = x*(1.-x + c1*x*x);
606 w = x*(1.0-x + c1*x*x);
607 xi = log(1.0-x + c1*x*x) - w;
614 w = -1 + sqrt(2.0*exp(1.))*sqrt(x-turnpt);
629 while ( relerr > mach_eps && i < maxit ) {
633 s = 6.*(w+1.0)*(w+1.0);
634 w = w * ( 1.0 + r + r2/(2.0*( w+1.0)) - (2. * w -1.0)*r3/s );
635 w = ((w*x < 0.0) ? -w : w);
638 relerr = ((x > 1.0) ? xi/w : xi);
639 relerr = ((relerr < 0.0) ? -relerr : relerr);
642 ierr = ((i == maxit) ? 2 : ierr);
658 Real arg1 = (Vsrc + Is*Rs)/
Vth_;
659 Real evd = exp(arg1);
660 Real lambWArg = Is*Rs*evd/
Vth_;
661 Real lambWReturn = 0.0;
662 Real lambWError = 0.0;
664 lambertw(lambWArg, lambWReturn, ierr, lambWError);
666 std::cout <<
"LambertW error: argument is not in the domain" << std::endl;
670 std::cout <<
"LambertW error: BUG!" << std::endl;
672 Real Id = -Is+
Vth_*(lambWReturn)/Rs;
687 ROL::Ptr<vector> lambdap =
getVector(lambda);
688 ROL::Ptr<const vector> Ip =
getVector(I);
689 ROL::Ptr<const vector> Sp =
getVector(S);
692 for (
uint i = 0; i < n; i++ ){
693 (*lambdap)[i] = ((*Imeas_)[i]-(*Ip)[i])
694 /
diodeI((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
708 ROL::Ptr<vector> sensp =
getVector(sens);
709 ROL::Ptr<const vector> Ip =
getVector(I);
710 ROL::Ptr<const vector> Sp =
getVector(S);
713 for (
uint i = 0; i < n; i++ ) {
714 (*sensp)[i] = -
diodeIs((*Ip)[i],(*
Vsrc_)[i],(*Sp)[0],(*Sp)[1])
729 ROL::Ptr<vector> sensp =
getVector(sens);
730 ROL::Ptr<const vector> Ip =
getVector(I);
731 ROL::Ptr<const vector> Sp =
getVector(S);
734 for (
uint i = 0; i < n; i++ ) {
735 (*sensp)[i] = -
diodeRs((*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.
typename PV< Real >::size_type size_type
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...
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...
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.
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
ROL::Ptr< const vector > getVector(const V &x)
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.
ROL::Ptr< vector > getVector(V &x)
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.
ROL::Ptr< std::vector< Real > > Vsrc_
Vector of source voltages in DC analysis (input)
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.
ROL::Ptr< std::vector< Real > > Imeas_
Vector of measured currents in DC analysis (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.
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.