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
18namespace ROL {
19namespace ZOO {
20
32template<class Real>
33class 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
42private:
44 Real Vth_;
46 ROL::Ptr<std::vector<Real> > Imeas_;
48 ROL::Ptr<std::vector<Real> > Vsrc_;
50 bool lambertw_;
52 Real noise_;
59
60public:
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
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
407private:
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 std::vector implementation of the ROL::Vector interface that handles scalings in the inn...
Provides the interface to evaluate objective functions.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Provides the std::vector implementation of the ROL::Vector interface that handles scalings in the inn...
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
The diode circuit problem.
void solve_sensitivity_Is(Vector< Real > &sens, const Vector< Real > &I, const Vector< Real > &S)
Solve the sensitivity equation wrt Is.
Real Vth_
Thermal voltage (constant)
Real diodeII(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Second derivative of diode equation wrt I^2.
Real value(const Vector< Real > &S, Real &tol)
Evaluate objective function.
Real diode(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Diode equation.
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 > > Vsrc_
Vector of source voltages in DC analysis (input)
Real diodeIIs(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Second derivative of diode equation wrt I and Is.
Real diodeIs(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Derivative of diode equation wrt Is.
Real diodeI(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Derivative of diode equation wrt I.
PrimalScaledStdVector< Real > PSV
DualScaledStdVector< Real > DSV
Real Newton(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Newton's method with line search.
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 random(const Real left, const Real right) const
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.
void solve_adjoint(Vector< Real > &lambda, const Vector< Real > &I, const Vector< Real > &S)
Solve the adjoint equation.
void set_method(bool lambertw)
Change the method for solving the circuit if needed.
Real lambertWCurrent(Real Is, Real Rs, Real Vsrc)
Find currents using Lambert-W function.
ROL::Ptr< const vector > getVector(const V &x)
Real diodeIRs(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Second derivative of diode equation wrt I and Rs.
Real noise_
Percentage of noise to add to measurements; if 0.0 - no noise.
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 diodeRs(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Derivative of diode equation wrt Rs
Real diodeIsRs(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Second derivative of diode equation wrt Is and Rs.
ROL::Ptr< std::vector< Real > > Imeas_
Vector of measured currents in DC analysis (data)
ROL::Ptr< vector > getVector(V &x)
bool lambertw_
If true, use Lambert-W function to solve circuit, else use Newton's method.
Real diodeRsRs(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Second derivative of diode equation wrt Rs^2.
void lambertw(Real x, Real &w, int &ierr, Real &xi)
Lambert-W function for diodes.
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.
bool use_adjoint_
If true, use adjoint gradient computation, else compute gradient using sensitivities.
void solve_sensitivity_Rs(Vector< Real > &sens, const Vector< Real > &I, const Vector< Real > &S)
Solve the sensitivity equation wrt Rs.
void solve_circuit(Vector< Real > &I, const Vector< Real > &S)
Solve circuit given optimization parameters Is and Rs.
Real diodeIsIs(const Real I, const Real Vsrc, const Real Is, const Real Rs)
Second derivative of diode equation wrt Is^2.
constexpr auto dim