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