ROL
ROL_DiodeCircuit.hpp
Go to the documentation of this file.
1 #ifndef ROL_DIODECIRCUIT_HPP
2 #define ROL_DIODECIRCUIT_HPP
3 
4 #include "ROL_Objective.hpp"
7 
8 #include <iostream>
9 #include <fstream>
10 #include <string>
11 
18 namespace ROL {
19 namespace ZOO {
20 
32 template<class Real>
33 class Objective_DiodeCircuit : public Objective<Real> {
34 
35  typedef std::vector<Real> vector;
36  typedef Vector<Real> V;
40  typedef typename vector::size_type uint;
41 
42 private:
44  Real Vth_;
46  ROL::Ptr<std::vector<Real> > Imeas_;
48  ROL::Ptr<std::vector<Real> > Vsrc_;
50  bool lambertw_;
52  Real noise_;
59 
60 public:
61 
76  Objective_DiodeCircuit(Real Vth, Real Vsrc_min, Real Vsrc_max, Real Vsrc_step,
77  Real true_Is, Real true_Rs,
78  bool lambertw, Real noise,
79  bool use_adjoint, int use_hessvec)
80  : Vth_(Vth), lambertw_(lambertw), use_adjoint_(use_adjoint), use_hessvec_(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;
86  // Generate problem data
87  if ( lambertw_ ) {
88  std::cout << "Generating data using Lambert-W function." << std::endl;
89  }
90  else {
91  std::cout << "Generating data using Newton's method." << std::endl;
92  }
93  for ( int i = 0; i < n; i++ ) {
94  (*Vsrc_)[i] = Vsrc_min+i*Vsrc_step;
95  if (lambertw_) {
96  (*Imeas_)[i] = lambertWCurrent(true_Is,true_Rs,(*Vsrc_)[i]);
97  }
98  else {
99  Real I0 = 1.e-12; // initial guess for Newton
100  (*Imeas_)[i] = Newton(I0,Vsrc_min+i*Vsrc_step,true_Is,true_Rs);
101  }
102  if ( noise > 0.0 ) {
103  (*Imeas_)[i] += noise*pow(10,(int)log10((*Imeas_)[i]))*random(left, right);
104  }
105  // Write generated data into file
106  if( output.is_open() ) {
107  output << std::setprecision(8) << std::scientific << (*Vsrc_)[i] << " " << (*Imeas_)[i] << "\n";
108  }
109  }
110  output.close();
111  }
112 
126  Objective_DiodeCircuit(Real Vth, std::ifstream& input_file,
127  bool lambertw, Real noise,
128  bool use_adjoint, int use_hessvec)
129  : Vth_(Vth), lambertw_(lambertw), use_adjoint_(use_adjoint), use_hessvec_(use_hessvec) {
130  std::string line;
131  int dim = 0;
132  for( int k = 0; std::getline(input_file,line); ++k ) {
133  dim = k;
134  } // count number of lines
135  input_file.clear(); // reset to beginning of file
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);
139  Real Vsrc, Imeas;
140  std::cout << "Using input file to generate data." << "\n";
141  for( int i = 0; i < dim; i++ ){
142  input_file >> Vsrc;
143  input_file >> Imeas;
144  (*Vsrc_)[i] = Vsrc;
145  (*Imeas_)[i] = Imeas;
146  }
147  input_file.close();
148  }
149 
151  void set_method(bool lambertw){
153  }
154 
157 
158  ROL::Ptr<vector> Ip = getVector(I);
159  ROL::Ptr<const vector> Sp = getVector(S);
160 
161  uint n = Ip->size();
162 
163  if ( lambertw_ ) {
164  // Using Lambert-W function
165  Real lambval;
166  for ( uint i = 0; i < n; i++ ) {
167  lambval = lambertWCurrent((*Sp)[0],(*Sp)[1],(*Vsrc_)[i]);
168  (*Ip)[i] = lambval;
169  }
170  }
171  else{
172  // Using Newton's method
173  Real I0 = 1.e-12; // Initial guess for Newton
174  for ( uint i = 0; i < n; i++ ) {
175  (*Ip)[i] = Newton(I0,(*Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
176  }
177  }
178  }
179 
187  Real value(const Vector<Real> &S, Real &tol){
188 
189  ROL::Ptr<const vector> Sp = getVector(S);
190  uint n = Imeas_->size();
191  STDV I( ROL::makePtr<vector>(n,0.0) );
192  ROL::Ptr<vector> Ip = getVector(I);
193 
194  // Solve state equation
195  solve_circuit(I,S);
196  Real val = 0;
197 
198  for ( uint i = 0; i < n; i++ ) {
199  val += ((*Ip)[i]-(*Imeas_)[i])*((*Ip)[i]-(*Imeas_)[i]);
200  }
201  return val/2.0;
202  }
203 
205  void gradient(Vector<Real> &g, const Vector<Real> &S, Real &tol){
206 
207 
208  ROL::Ptr<vector> gp = getVector(g);
209  ROL::Ptr<const vector> Sp = getVector(S);
210 
211  uint n = Imeas_->size();
212 
213  STDV I( ROL::makePtr<vector>(n,0.0) );
214  ROL::Ptr<vector> Ip = getVector(I);
215 
216  // Solve state equation
217  solve_circuit(I,S);
218 
219  if ( use_adjoint_ ) {
220  // Compute the gradient of the reduced objective function using adjoint computation
221  STDV lambda( ROL::makePtr<vector>(n,0.0) );
222  ROL::Ptr<vector> lambdap = getVector(lambda);
223 
224  // Solve adjoint equation
225  solve_adjoint(lambda,I,S);
226 
227  // Compute gradient
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];
232  }
233  }
234  else {
235  // Compute the gradient of the reduced objective function using sensitivities
236  STDV sensIs( ROL::makePtr<vector>(n,0.0) );
237  STDV sensRs( ROL::makePtr<vector>(n,0.0) );
238  // Solve sensitivity equations
239  solve_sensitivity_Is(sensIs,I,S);
240  solve_sensitivity_Rs(sensRs,I,S);
241 
242  ROL::Ptr<vector> sensIsp = getVector(sensIs);
243  ROL::Ptr<vector> sensRsp = getVector(sensRs);
244 
245  // Write sensitivities into file
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";
250  }
251  }
252  output.close();
253  // Compute gradient
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];
258  }
259  }
260  }
261 
269  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &S, Real &tol ){
270 
271 
272 
273  if ( use_hessvec_ == 0 ) {
274  Objective<Real>::hessVec(hv, v, S, tol);
275  }
276  else if ( use_hessvec_ == 1 ) {
277  ROL::Ptr<vector> hvp = getVector(hv);
278  ROL::Ptr<const vector> vp = getVector(v);
279  ROL::Ptr<const vector> Sp = getVector(S);
280 
281  uint n = Imeas_->size();
282 
283  STDV I( ROL::makePtr<vector>(n,0.0) );
284  ROL::Ptr<vector> Ip = getVector(I);
285 
286  // Solve state equation
287  solve_circuit(I,S);
288 
289  STDV lambda( ROL::makePtr<vector>(n,0.0) );
290  ROL::Ptr<vector> lambdap = getVector(lambda);
291 
292  // Solve adjoint equation
293  solve_adjoint(lambda,I,S);
294 
295  STDV w( ROL::makePtr<vector>(n,0.0) );
296  ROL::Ptr<vector> wp = getVector(w);
297 
298  // Solve state sensitivity equation
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]);
303  }
304 
305  STDV p( ROL::makePtr<vector>(n,0.0) );
306  ROL::Ptr<vector> pp = getVector(p);
307 
308  // Solve for p
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] );
314  }
315 
316  // Assemble Hessian-vector product
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] );
327  }
328  }
329  else if ( use_hessvec_ == 2 ) {
330  //Gauss-Newton approximation
331  ROL::Ptr<vector> hvp = getVector(hv);
332  ROL::Ptr<const vector> vp = getVector(v);
333  ROL::Ptr<const vector> Sp = getVector(S);
334 
335  uint n = Imeas_->size();
336 
337  STDV I( ROL::makePtr<vector>(n,0.0) );
338  ROL::Ptr<vector> Ip = getVector(I);
339 
340  // Solve state equation
341  solve_circuit(I,S);
342 
343  // Compute sensitivities
344  STDV sensIs( ROL::makePtr<vector>(n,0.0) );
345  STDV sensRs( ROL::makePtr<vector>(n,0.0) );
346 
347  // Solve sensitivity equations
348  solve_sensitivity_Is(sensIs,I,S);
349  solve_sensitivity_Rs(sensRs,I,S);
350  ROL::Ptr<vector> sensIsp = getVector(sensIs);
351  ROL::Ptr<vector> sensRsp = getVector(sensRs);
352 
353  // Compute approximate Hessian
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];
359  }
360 
361  // Compute approximate Hessian-times-vector
362  (*hvp)[0] = H11*(*vp)[0] + H12*(*vp)[1];
363  (*hvp)[1] = H12*(*vp)[0] + H22*(*vp)[1];
364  }
365  else {
366  ROL::Objective<Real>::hessVec( hv, v, S, tol ); // Use parent class function
367  }
368  }
369 
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);
383  StdVector<Real> S(S_ptr);
384  std::ofstream output ("Objective.dat");
385 
386  Real Is = 0.0;
387  Real Rs = 0.0;
388  Real val = 0.0;
389  Real tol = 1.e-16;
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;
396  (*S_ptr)[0] = Is;
397  (*S_ptr)[1] = Rs;
398  val = value(S,tol);
399  if ( output.is_open() ) {
400  output << std::scientific << Is << " " << Rs << " " << val << std::endl;
401  }
402  }
403  }
404  output.close();
405  }
406 
407 private:
408 
409  ROL::Ptr<const vector> getVector( const V& x ) {
410  try {
411  return dynamic_cast<const STDV&>(x).getVector();
412  }
413  catch (std::exception &e) {
414  try {
415  return dynamic_cast<const PSV&>(x).getVector();
416  }
417  catch (std::exception &e) {
418  return dynamic_cast<const DSV&>(x).getVector();
419  }
420  }
421  }
422 
423  ROL::Ptr<vector> getVector( V& x ) {
424 
425  try {
426  return dynamic_cast<STDV&>(x).getVector();
427  }
428  catch (std::exception &e) {
429  try {
430  return dynamic_cast<PSV&>(x).getVector();
431  }
432  catch (std::exception &e) {
433  return dynamic_cast<DSV&>(x).getVector();
434  }
435  }
436  }
437 
438  Real random(const Real left, const Real right) const {
439  return (Real)rand()/(Real)RAND_MAX * (right - left) + left;
440  }
441 
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);
454  }
455 
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_);
459  }
460 
462  Real diodeIs(const Real I, const Real Vsrc, const Real Is, const Real Rs){
463  return 1-exp((Vsrc-I*Rs)/Vth_);
464  }
465 
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_);
469  }
470 
472  Real diodeII(const Real I, const Real Vsrc, const Real Is, const Real Rs){
473  return -Is*exp((Vsrc-I*Rs)/Vth_)*(Rs/Vth_)*(Rs/Vth_);
474  }
475 
477  Real diodeIIs(const Real I, const Real Vsrc, const Real Is, const Real Rs){
478  return exp((Vsrc-I*Rs)/Vth_)*(Rs/Vth_);
479  }
480 
482  Real diodeIRs(const Real I, const Real Vsrc, const Real Is, const Real Rs){
483  return (Is/Vth_)*exp((Vsrc-I*Rs)/Vth_)*(1-(I*Rs)/Vth_);
484  }
485 
487  Real diodeIsIs(const Real I, const Real Vsrc, const Real Is, const Real Rs){
488  return 0;
489  }
490 
492  Real diodeIsRs(const Real I, const Real Vsrc, const Real Is, const Real Rs){
493  return exp((Vsrc-I*Rs)/Vth_)*(I/Vth_);
494  }
495 
497  Real diodeRsRs(const Real I, const Real Vsrc, const Real Is, const Real Rs){
498  return -Is*exp((Vsrc-I*Rs)/Vth_)*(I/Vth_)*(I/Vth_);
499  }
500 
508  Real Newton(const Real I, const Real Vsrc, const Real Is, const Real Rs){
509  Real EPS = 1.e-16;
510  Real TOL = 1.e-13;
511  int MAXIT = 200;
512  Real IN = I;
513  Real fval = diode(IN,Vsrc,Is,Rs);
514  Real dfval = 0.0;
515  Real IN_tmp = 0.0;
516  Real fval_tmp = 0.0;
517  Real alpha = 1.0;
518  for ( int i = 0; i < MAXIT; i++ ) {
519  if ( std::abs(fval) < TOL ) {
520  // std::cout << "converged with |fval| = " << std::abs(fval) << " and TOL = " << TOL << "\n";
521  break;
522  }
523  dfval = diodeI(IN,Vsrc,Is,Rs);
524  if( std::abs(dfval) < EPS ){
525  std::cout << "denominator is too small" << std::endl;
526  break;
527  }
528 
529  alpha = 1.0;
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) ) {
533  alpha /= 2.0;
534  IN_tmp = IN - alpha*fval/dfval;
535  fval_tmp = diode(IN_tmp,Vsrc,Is,Rs);
536  if ( alpha < std::sqrt(EPS) ) {
537  // std::cout << "Step tolerance met\n";
538  break;
539  }
540  }
541  IN = IN_tmp;
542  fval = fval_tmp;
543  // if ( i == MAXIT-1){
544  // std::cout << "did not converge " << std::abs(fval) << "\n";
545  // }
546  }
547  return IN;
548  }
549 
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;
585  mach_eps = 2.e-15; // float:2e-7
586  ierr = 0;
587 
588  if ( x > c1 ) {
589  w = c2*log(x);
590  xi = log( x/ w) - w;
591  }
592  else {
593  if ( x >= 0.0 ) {
594  w = x;
595  if ( x == 0. ) {
596  return;
597  }
598  if ( x < (1-c2) ) {
599  w = x*(1.-x + c1*x*x);
600  }
601  xi = - w;
602  }
603  else {
604  if ( x >= turnpt ){
605  if ( x > -0.2 ){
606  w = x*(1.0-x + c1*x*x);
607  xi = log(1.0-x + c1*x*x) - w;
608  }
609  else {
610  diff = x-turnpt;
611  if ( diff < 0.0 ) {
612  diff = -diff;
613  }
614  w = -1 + sqrt(2.0*exp(1.))*sqrt(x-turnpt);
615  if ( diff == 0.0 ) {
616  return;
617  }
618  xi = log( x/ w) - w;
619  }
620  }
621  else {
622  ierr = 1; // x is not in the domain.
623  w = -1.0;
624  return;
625  }
626  }
627  }
628 
629  while ( relerr > mach_eps && i < maxit ) {
630  r = xi/(w+1.0); //singularity at w=-1
631  r2 = r*r;
632  r3 = r2*r;
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);
636  xi = log( x/ w) - w;
637 
638  relerr = ((x > 1.0) ? xi/w : xi);
639  relerr = ((relerr < 0.0) ? -relerr : relerr);
640  ++i;
641  }
642  ierr = ((i == maxit) ? 2 : ierr);
643  }
644 
657  Real lambertWCurrent(Real Is, Real Rs, Real Vsrc){
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;
663  int ierr = 0;
664  lambertw(lambWArg, lambWReturn, ierr, lambWError);
665  if ( ierr == 1 ) {
666  std::cout << "LambertW error: argument is not in the domain" << std::endl;
667  return -1.0;
668  }
669  if ( ierr == 2 ) {
670  std::cout << "LambertW error: BUG!" << std::endl;
671  }
672  Real Id = -Is+Vth_*(lambWReturn)/Rs;
673  //Real Gd = lambWReturn / ((1 + lambWReturn)*RS);
674  return Id;
675  }
676 
684  void solve_adjoint(Vector<Real> &lambda, const Vector<Real> &I, const Vector<Real> &S){
685 
686 
687  ROL::Ptr<vector> lambdap = getVector(lambda);
688  ROL::Ptr<const vector> Ip = getVector(I);
689  ROL::Ptr<const vector> Sp = getVector(S);
690 
691  uint n = Ip->size();
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]);
695  }
696  }
697 
706 
707 
708  ROL::Ptr<vector> sensp = getVector(sens);
709  ROL::Ptr<const vector> Ip = getVector(I);
710  ROL::Ptr<const vector> Sp = getVector(S);
711 
712  uint n = Ip->size();
713  for ( uint i = 0; i < n; i++ ) {
714  (*sensp)[i] = -diodeIs((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1])
715  /diodeI((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
716  }
717  }
718 
727 
728 
729  ROL::Ptr<vector> sensp = getVector(sens);
730  ROL::Ptr<const vector> Ip = getVector(I);
731  ROL::Ptr<const vector> Sp = getVector(S);
732 
733  uint n = Ip->size();
734  for ( uint i = 0; i < n; i++ ) {
735  (*sensp)[i] = -diodeRs((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1])
736  /diodeI((*Ip)[i],(*Vsrc_)[i],(*Sp)[0],(*Sp)[1]);
737  }
738  }
739 }; // class Objective_DiodeCircuit
740 
741 
742  // template<class Real>
743  // void getDiodeCircuit( ROL::Ptr<Objective<Real> > &obj, Vector<Real> &x0, Vector<Real> &x ) {
744  // // Cast Initial Guess and Solution Vectors
745  // ROL::Ptr<std::vector<Real> > x0p =
746  // ROL::constPtrCast<std::vector<Real> >((dynamic_cast<PrimalScaledStdVector<Real>&>(x0)).getVector());
747  // ROL::Ptr<std::vector<Real> > xp =
748  // ROL::constPtrCast<std::vector<Real> >((dynamic_cast<PrimalScaledStdVector<Real>&>(x)).getVector());
749 
750  // int n = xp->size();
751 
752  // // Resize Vectors
753  // n = 2;
754  // x0p->resize(n);
755  // xp->resize(n);
756 
757  // // Instantiate Objective Function
758  // obj = ROL::makePtr<Objective_DiodeCircuit<Real>>(0.02585,0.0,1.0,1.e-2);
759  // //ROL::Objective_DiodeCircuit<Real> obj(0.02585,0.0,1.0,1.e-2);
760 
761  // // Get Initial Guess
762  // (*x0p)[0] = 1.e-13;
763  // (*x0p)[1] = 0.2;
764 
765  // // Get Solution
766  // (*xp)[0] = 1.e-12;
767  // (*xp)[1] = 0.25;
768 
769  // }
770 
771 
772 } //end namespace ZOO
773 } //end namespace ROL
774 
775 #endif
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.
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.
Definition: ROL_Vector.hpp:46
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&#39;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&#39;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.
constexpr auto dim
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.