diff --git a/.gitignore b/.gitignore index f18eba0007..cbbc8257e6 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ # Byte-compiled / optimized / DLL files +.pytest_cache/ __pycache__/ *.py[cod] diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst index 5909db3a7a..31cfabdb02 100644 --- a/docs/sphinx/source/api.rst +++ b/docs/sphinx/source/api.rst @@ -65,7 +65,6 @@ algorithm. spa - Correlations and analytical expressions for low precision solar position calculations. @@ -80,6 +79,7 @@ calculations. solarposition.equation_of_time_pvcdrom solarposition.hour_angle + Clear sky ========= @@ -206,6 +206,18 @@ Functions relevant for the single diode model. pvsystem.i_from_v pvsystem.singlediode pvsystem.v_from_i + pvsystem.max_power_point + +Low-level functions for solving the single diode equation. + +.. autosummary:: + :toctree: generated/ + + singlediode_methods.estimate_voc + singlediode_methods.bishop88 + singlediode_methods.bishop88_i_from_v + singlediode_methods.bishop88_v_from_i + singlediode_methods.bishop88_mpp SAPM model ---------- @@ -232,7 +244,6 @@ PVWatts model pvsystem.pvwatts_ac pvsystem.pvwatts_losses - Other ----- diff --git a/docs/sphinx/source/whatsnew/v0.6.0.rst b/docs/sphinx/source/whatsnew/v0.6.0.rst index 5f8d938a98..009f015e62 100644 --- a/docs/sphinx/source/whatsnew/v0.6.0.rst +++ b/docs/sphinx/source/whatsnew/v0.6.0.rst @@ -11,10 +11,43 @@ API Changes Enhancements ~~~~~~~~~~~~ -* Add sea surface albedo in irradiance.py (:issue:`458`) -* Implement first_solar_spectral_loss in modelchain.py (:issue:`359`) -* Clarify arguments Egref and dEgdT for calcparams_desoto (:issue:`462`) -* Add pvsystem.calcparams_pvsyst to compute values for the single diode equation using the PVsyst v6 model (:issue:'470') +* Add sea surface albedo in ``irradiance.py`` (:issue:`458`) +* Implement :meth:`~pvlib.modelchain.ModelChain.first_solar_spectral_loss` + in ``modelchain.py`` (:issue:`359`) +* Clarify arguments ``Egref`` and ``dEgdT`` for + :func:`~pvlib.pvsystem.calcparams_desoto` (:issue:`462`) +* Add pvsystem.calcparams_pvsyst to compute values for the single diode equation + using the PVsyst v6 model (:issue:'470') +* Extend :func:`~pvlib.pvsystem.singlediode` with an additional keyword argument + ``method`` in ``('lambertw', 'newton', 'brentq')``, default is ``'lambertw'``, + to select a method to solve the single diode equation for points on the IV + curve. Selecting either ``'brentq'`` or ``'newton'`` as the method uses + :func:`~pvlib.singlediode_methods.bishop88` with the corresponding method. + (:issue:`410`) +* Implement new methods ``'brentq'`` and ``'newton'`` for solving the single + diode equation for points on the IV curve. ``'brentq'`` uses a bisection + method (Brent, 1973) that may be slow but guarantees a solution. ``'newton'`` + uses the Newton-Raphson method and may be faster but is not guaranteed to + converge. However, ``'newton'`` should be safe for well-behaved IV curves. + (:issue:`408`) +* Implement :func:`~pvlib.singlediode_methods.bishop88` for explicit calculation + of arbitrary IV curve points using diode voltage instead of cell voltage. If + ``method`` is either ``'newton'`` or ``'brentq'`` and ``ivcurve_pnts`` in + :func:`~pvlib.pvsystem.singlediode` is provided, the IV curve points will be + log spaced instead of linear. +* Implement :func:`~pvlib.singlediode_methods.estimate_voc` to estimate open + circuit voltage by assuming :math:`R_{sh} \to \infty` and :math:`R_s=0` as an + upper bound in bisection method for :func:`~pvlib.pvsystem.singlediode` when + method is either ``'newton'`` or ``'brentq'``. +* Add :func:`~pvlib.pvsystem.max_power_point` method to compute the max power + point using the new ``'brentq'`` method. +* Add new module ``pvlib.singlediode_methods`` with low-level functions for + solving the single diode equation such as: + :func:`~pvlib.singlediode_methods.bishop88`, + :func:`~pvlib.singlediode_methods.estimate_voc`, + :func:`~pvlib.singlediode_methods.bishop88_i_from_v`, + :func:`~pvlib.singlediode_methods.bishop88_v_from_i`, and + :func:`~pvlib.singlediode_methods.bishop88_mpp`. Bug fixes @@ -50,5 +83,6 @@ Contributors * Will Holmgren * Yu Cao * Cliff Hansen +* Mark Mikofski * Alan Mathew diff --git a/pvlib/__init__.py b/pvlib/__init__.py index b756147eba..3a231da728 100644 --- a/pvlib/__init__.py +++ b/pvlib/__init__.py @@ -11,3 +11,4 @@ from pvlib import pvsystem from pvlib import spa from pvlib import modelchain +from pvlib import singlediode_methods diff --git a/pvlib/data/bishop88_numerical_precision.csv b/pvlib/data/bishop88_numerical_precision.csv new file mode 100644 index 0000000000..94cf44adc9 --- /dev/null +++ b/pvlib/data/bishop88_numerical_precision.csv @@ -0,0 +1,101 @@ +,grad,grad2p,grad_i,grad_p,grad_v,i,p,v +0.0,-0.0029752937017359805,-0.005957140772172728,-0.0029785726893041607,5.870505570369785,1.0011020718950425,5.86405008,-12.723220836076763,-2.1696985296 +0.5845170900989292,-0.00297529507385398,-0.005957144029576585,-0.002978574064448173,5.8670235188728626,1.0011020724038457,5.862309052970334,-9.289047121136205,-1.5845372595000944 +1.1690341801978583,-0.002975296745849822,-0.005957148329744609,-0.0029785757401313675,5.863541465180093,1.0011020730238487,5.860568025051929,-5.856910966556231,-0.9993759890713553 +1.7535512702967875,-0.002975298783262126,-0.0059571539728827455,-0.002978577782036902,5.860059408597745,1.0011020737793537,5.85882699605055,-2.426812373397211,-0.41421471824191614 +2.3380683603957166,-0.002975301265953046,-0.005957161340610417,-0.0029785802702030567,5.856577348233737,1.001102074699975,5.857085965729514,1.0012486569100834,0.17094655307579656 +2.922585450494646,-0.0029753042912386935,-0.005957170917207319,-0.002978583302160564,5.8530952829440785,1.0011020758217994,5.855344933800406,4.427272122453232,0.7561078249884956 +3.507102540593575,-0.0029753079777037196,-0.005957183316258555,-0.0029785869967556004,5.849613211265367,1.0011020771887995,5.853603899911782,7.851258020695973,1.3412690976262158 +4.091619630692504,-0.0029753124698495737,-0.005957199314042896,-0.0029785914988083007,5.846131131329814,1.001102078854559,5.851862863635386,11.27320634829687,1.9264303711474113 +4.676136720791433,-0.002975317943758649,-0.005957219891338891,-0.002978596984789394,5.84264904075839,1.001102080884372,5.85012182444937,14.693117100878508,2.5115916457451664 +5.260653810890362,-0.0029753246139963283,-0.0059572462857335464,-0.002978603669737475,5.839166936526567,1.0011020833578028,5.8483807817178315,18.110990272734274,3.096752921654765 +5.845170900989292,-0.002975332742021483,-0.005957280057027362,-0.002978611815688049,5.835684814795811,1.0011020863718045,5.846639734665894,21.52682585645643,3.681914199162911 +6.429687991088221,-0.0029753426464350914,-0.005957323168961724,-0.002978621741944762,5.8322026707022365,1.0011020900445196,5.844898682349329,24.940623842465143,4.267075478618969 +7.01420508118715,-0.002975354715468701,-0.005957378091279118,-0.0029786338375954085,5.8287204980917755,1.0011020945199103,5.843157623617547,28.352384218412883,4.852236760448657 +7.598722171286079,-0.0029753694222022617,-0.005957447927099995,-0.0029786485767633332,5.825238289188552,1.0011020999734024,5.841416557068495,31.762106968432207,5.4373980451707355 +8.183239261385008,-0.0029753873431078255,-0.0059575365718074124,-0.002978666537192042,5.821756034179936,1.001102106618761,5.8396754809937015,35.16979207218701,6.022559333417338 +8.767756351483937,-0.002975409180645987,-0.005957648911127437,-0.0029786884228914967,5.818273720697683,1.0011021147164698,5.837934393311305,38.57543950367713,6.607720625958756 +9.352273441582867,-0.0029754357908008,-0.005957791067948803,-0.002978715091733781,5.81479133316957,1.0011021245839415,5.836193291484451,41.979049229733924,7.192881923733619 +9.936790531681796,-0.0029754682166324635,-0.0059579707097246176,-0.0029787475890798343,5.811308852009686,1.0011021366079595,5.834452172421854,45.38062120812886,7.778043227885708 +10.521307621780725,-0.0029755077291629788,-0.005958197431147377,-0.002978787188755338,5.8078262526078515,1.0011021512598395,5.832711032356624,48.78015538519764,8.363204539808773 +11.105824711879654,-0.0029755558771973863,-0.005958483230316493,-0.0029788354429819267,5.804343504069003,1.0011021691139033,5.830969866698612,52.177651692858575,8.948365861201166 +11.690341801978583,-0.0029756145480334634,-0.005958843100985721,-0.00297889424322092,5.800860567641557,1.0011021908699917,5.829228669854491,55.57311004487415,9.533527194132422 +12.274858892077512,-0.0029756860414395505,-0.005959295768885225,-0.0029789658943145158,5.797377394759002,1.0011022173808963,5.827487435008501,58.966530332166734,10.118688541124365 +12.859375982176442,-0.0029757731598002535,-0.00595986460680504,-0.002979053204830618,5.793893924600721,1.0011022496857873,5.825746153855294,62.35791241695444,10.703849905249983 +13.44389307227537,-0.0029758793179634957,-0.005960578771406571,-0.002979159597152625,5.790410081055457,1.0011022890509464,5.82400481627438,65.74725612541428,11.289011290253852 +14.0284101623743,-0.002976008677094632,-0.005961474614971044,-0.002979289241629458,5.786925768942823,1.001102337019403,5.822263409933423,69.13456123850973,11.874172700698933 +14.612927252473229,-0.0029761663077843346,-0.005962597437960581,-0.002979447220044218,5.783440869313562,1.0011023954714164,5.820521919804842,72.5198274805304,12.459334142145439 +15.197444342572158,-0.002976358388803591,-0.005964003663929506,-0.0029796397248090935,5.779955233606391,1.0011024666981794,5.818780327576764,75.90305450478175,13.044495621368755 +15.781961432671087,-0.002976592449296402,-0.005965763537684524,-0.002979874301694524,5.776468676386092,1.001102553491627,5.817038610935235,79.28424187572728,13.629657146625048 +16.366478522770016,-0.002976877663903328,-0.0059679644715214695,-0.0029801601456070836,5.7729809663218195,1.0011026592538745,5.815296742689563,82.66338904671491,14.214818727974878 +16.950995612868944,-0.0029772252123836693,-0.005970715193935012,-0.0029805084610099305,5.769491814983297,1.0011027881305736,5.8135546897065025,86.04049533221146,14.799980377677539 +17.535512702967875,-0.002977648717832082,-0.005974150891728205,-0.0029809329011135314,5.766000862931994,1.001102945173412,5.8118124116115,89.4155598732081,15.385142110671623 +18.120029793066806,-0.00297816478066586,-0.005978439581571335,-0.002981450103051982,5.762507662460098,1.0011031365381293,5.810069859206104,92.78858159413763,15.970303945160547 +18.704546883165733,-0.0029787936293126675,-0.005983790002784571,-0.002982080340022682,5.759011656176353,1.0011033697258085,5.808326972539485,96.15955914924689,16.555465903326123 +19.28906397326466,-0.0029795599131022697,-0.005990461391925088,-0.0029828483159518274,5.755512150447955,1.0011036538769023,5.806583678558498,99.5284908558716,17.140628012198015 +19.87358106336359,-0.0029804936684389004,-0.005998775584701608,-0.0029837841338348874,5.752008282472884,1.0011040001295188,5.804839888244158,102.89537461145103,17.72579030471325 +20.458098153462522,-0.002981631496121704,-0.006009131995583255,-0.0029849244757089214,5.748498979467168,1.0011044220560124,5.8030954931222825,106.26020779036207,18.310952821007277 +21.04261524356145,-0.002983017995955065,-0.006022026154854778,-0.0029863140405090977,5.744982908093305,1.0011049361949884,5.8013503610115436,109.62298711571498,18.896115609987177 +21.627132333660377,-0.0029847075148726693,-0.006038072642512993,-0.0029880072961702775,5.741458411813756,1.001105562699583,5.799604330842244,112.98370850009512,19.481278731248747 +22.211649423759308,-0.0029867662770834236,-0.006058033455340477,-0.002990070614652262,5.73792343330725,1.0011063261274213,5.797857206342731,116.34236684779933,20.066442257412493 +22.79616651385824,-0.002989274979714781,-0.00608285308640337,-0.002992584873577066,5.734375418411309,1.0011072564032235,5.7961087483459774,119.69895580934383,20.651606276970227 +23.380683603957166,-0.0029923319556650167,-0.00611370189577642,-0.0029956486264567862,5.730811197222261,1.001108389991789,5.794358665414746,123.05346747682691,21.23677089775371 +23.965200694056094,-0.00299605702759383,-0.006152029720641873,-0.0029993819657781226,5.727226836956887,1.0011097713273378,5.792606602417884,126.40589200602012,21.82193625116148 +24.549717784155025,-0.003000596204048713,-0.0061996321282111674,-0.003003931230368128,5.723617459912789,1.0011114545552362,5.790852126609946,129.75621714771086,22.40710249730935 +25.134234874253956,-0.0030061274017001574,-0.00625873227607497,-0.0030094747415597284,5.719977018301338,1.0011135056543772,5.789094710668513,133.1044276666803,22.992269831306608 +25.718751964352883,-0.0030128674178291435,-0.006332082036065065,-0.003016229793002312,5.716298015799175,1.001116005023411,5.7873337120242905,136.45050462158724,23.577438490903894 +26.30326905445181,-0.003021080426142027,-0.006423086889588923,-0.0030244611681028225,5.712571163286861,1.001119050632198,5.785568347673785,139.79442447271,24.16260876581251 +26.88778614455074,-0.003031088328586912,-0.006535960151633435,-0.0030344915189626518,5.708784953312269,1.0011227618620162,5.783797663487256,143.13615797669814,24.747781009060454 +27.472303234649672,-0.0030432833684303404,-0.006675913372617289,-0.0030467140136421905,5.704925134203544,1.0011272841850476,5.782020496808899,146.4756688178521,25.332955650830378 +28.0568203247486,-0.003058143498242408,-0.006849391357709309,-0.0030616077474984646,5.700974060304231,1.0011327948665745,5.780235430883248,149.81291191355504,25.9181332153218 +28.641337414847527,-0.003076251104057797,-0.007064362200442791,-0.003079756522686961,5.696909889317819,1.0011395099133942,5.77844073932143,153.14783131680173,26.503314341298598 +29.225854504946458,-0.003098315817990082,-0.007330675135155671,-0.0031018717319434906,5.692705590992404,1.0011476925408191,5.776634318430448,156.4803576206475,27.088499807127196 +29.81037159504539,-0.003125202311026061,-0.007660501973435053,-0.003128820243640859,5.688327723056219,1.0011576634901471,5.774813604752972,159.81040474704287,27.67369056128679 +30.394888685144316,-0.0031579641517379847,-0.008068881528586546,-0.0031616583811542435,5.683734920072716,1.001169813601027,5.772975474585356,163.13786597492938,28.258887759547736 +30.979405775243244,-0.0031978850526391375,-0.00857439090203128,-0.00320167332845292,5.678876028280287,1.0011846191315277,5.7711161215352025,166.46260902843693,28.844092810275217 +31.563922865342175,-0.003246529112840734,-0.009199972992622521,-0.003250433584927979,5.673687803977977,1.0012026604264233,5.769230907319024,169.78447000403867,29.42930742963414 +32.148439955441106,-0.0033058020143741994,-0.009973956320365261,-0.0033098504471761533,5.668092073956306,1.0012246446654551,5.767314179951555,173.10324586373662,30.01453370885903 +32.73295704554003,-0.003378025553070552,-0.010931311503920758,-0.0033822529276927283,5.661992233047112,1.0012514335832463,5.765359052200152,176.41868515749206,30.599774196225976 +33.31747413563896,-0.0034660283989182195,-0.012115198827249804,-0.0034704790471242375,5.655268925097611,1.001284077247436,5.763357131620177,179.73047655936463,31.185031996939493 +33.90199122573789,-0.0035732566041354074,-0.013578873671946769,-0.003577987078539832,5.647774718368853,1.0013238552190598,5.761298191589311,183.03823470472858,31.770310894849842 +34.48650831583682,-0.0037039081322002257,-0.015388031652369478,-0.0037089911042561515,5.6393275430704275,1.0013723267085748,5.75916977044608,186.34148269622136,32.35561550077176 +35.07102540593575,-0.0038630965943735403,-0.017623693632213337,-0.00386862619875006,5.629702605722588,1.0014313916935376,5.756956683019679,189.63963049849377,32.94095143321847 +35.65554249603468,-0.004057050483175599,-0.02038575308226625,-0.004063149712472153,5.61862243017849,1.0015033653936147,5.754640425404113,192.93194825988823,33.52632553863516 +36.24005958613361,-0.004293355521479788,-0.023797335221513505,-0.00430018654645125,5.605744595922764,1.001591069022187,5.752198449645108,196.21753337488724,34.11174615976492 +36.82457667623254,-0.004581249343709803,-0.02801014993255118,-0.004589028031923824,5.5906466476921715,1.0016979403718118,5.7496032799090635,199.4952698246832,34.69722346266619 +37.409093766331466,-0.004931979640265873,-0.033211059514087614,-0.004940996130423489,5.572807533010409,1.0018281685682566,5.746821435489804,202.7637779923818,35.28276983520024 +37.993610856430394,-0.005359239181470032,-0.039630128945234155,-0.005369887230186346,5.551584781783767,1.001986858275169,5.7438121184373045,206.02135272915035,35.86840037260859 +38.57812794652932,-0.005879693851065687,-0.04755048150237668,-0.005892512934731463,5.526186469961425,1.0021802297858506,5.740525614366289,209.26588692954923,36.454133469213794 +39.162645036628255,-0.006513623021513256,-0.057320347186224224,-0.00652935904135924,5.495636802151911,1.0024158628453028,5.736901343759845,212.49477723545002,37.0399915394371 +39.74716212672718,-0.007285695351029681,-0.06936776606936956,-0.007305388540099587,5.458733900288981,1.0027029937598368,5.732865487383382,215.7048077000665,37.62600189639534 +40.33167921682611,-0.008225907421738,-0.08421849332264997,-0.008251020108931858,5.413998091061136,1.0030528774403047,5.72832809273029,218.89200627168503,38.2121978225159 +40.916196306925045,-0.009370717591804082,-0.10251774615968948,-0.009403320460173364,5.359608638285938,1.003479228570264,5.723179548078225,222.05146775731706,38.7986198741361 +41.50071339702397,-0.010764412980085559,-0.12505653220817234,-0.010807457275435202,5.293326463334625,1.003998759191911,5.717286285946656,225.17713544588256,39.38531747122371 +42.0852304871229,-0.012460753541850314,-0.15280339779969054,-0.01251846968105451,5.212399935303023,1.0046318337819902,5.710485547540512,228.2615317418191,39.97235083453291 +42.66974757722183,-0.014524943507684662,-0.18694252257421318,-0.014603425662809928,5.113450295921546,1.0054032674952396,5.702579002957414,231.29542589991897,40.559793346127584 +43.254264667320754,-0.01703598663125314,-0.22891914566857377,-0.0171440509859117,4.992332723500172,1.0063432988647874,5.693324977084448,234.2674241568419,41.147734425799506 +43.83878175741969,-0.02008948701834557,-0.28049331005226186,-0.020239932668244195,4.84396846020808,1.0074887750872503,5.682428976456426,237.16346409409655,41.73628303613082 +44.423298847518616,-0.02380096064790911,-0.34380281144062375,-0.02401242257608301,4.662142873281297,1.0088845963531508,5.669532145748974,239.96619077809814,42.32557195359149 +45.00781593761754,-0.02830972227370975,-0.4214359713240549,-0.028609394154798077,4.439263870493397,1.0105854758372752,5.65419720142593,242.65418689950377,42.915762973089954 +45.59233302771648,-0.03378340555144074,-0.5165143255015552,-0.034211038748095036,4.166074866878158,1.0126580843367952,5.635891291170475,245.2010225160716,43.5070532499834 +46.176850117815405,-0.04042315708918318,-0.6327843977180755,-0.04103692870897761,3.8313166925969986,1.0151836636223217,5.613965107226749,247.57408175974606,44.099683028141506 +46.76136720791433,-0.04846951217513105,-0.7747162340807431,-0.04935462416110098,3.4213337208458725,1.0182612109396074,5.587627434940117,249.73311357612394,44.69394505698649 +47.34588429801326,-0.05820890367283452,-0.9476040789303638,-0.05949016077694885,2.9196214786872496,1.0220113594874711,5.555914138853929,251.6284406780968,45.29019606663731 +47.93040138811219,-0.06998066609906431,-1.1576612037329261,-0.07184082967105557,2.306316628998847,1.0265811069782906,5.517650370684814,253.19874470933536,45.88887075095881 +48.51491847821112,-0.08418426188957046,-1.4120961692589842,-0.0868907503524476,1.5576361908975433,1.0321495776304057,5.471404517810895,254.36832520582863,46.49049880662108 +49.09943556831005,-0.1012862621448542,-1.7191514760865074,-0.10522984716235959,0.6452820683996454,1.038935043450073,5.415432087153262,255.04370410222978,47.09572569606333 +49.68395265840898,-0.12182634556357387,-2.088077606107995,-0.12757697303254226,-0.46415966128739167,1.0472034800220407,5.347607324820204,255.10941464495394,47.705337948225505 +50.26846974850791,-0.14642122666854582,-2.5290062924988645,-0.15480808696547316,-1.8097768276699933,1.057278992177225,5.265339891152171,254.42277149133224,48.32029398878161 +50.85298683860684,-0.17576498907130353,-3.052677727102656,-0.18799058973276161,-3.436774886789651,1.0695565182011217,5.165473325009417,252.8073645834666,48.941761708353354 +51.437503928705766,-0.21062380604954392,-3.669969935228677,-0.22842516367685428,-5.396705965486928,1.084517310560436,5.044161317322018,250.04494912056487,49.57116424129662 +52.02202101880469,-0.25182254424459144,-4.391179275391126,-0.2776967566457863,-7.74735236946439,1.102747799958941,4.8967169440929945,245.8653121614329,50.21023574949029 +52.60653810890362,-0.3002203929971172,-5.225015825085368,-0.33773670852095095,-10.552107298134041,1.1249625821527518,4.717428949114363,239.93357550835128,50.861089397731305 +53.191055199002555,-0.3566726468607908,-6.177315083784205,-0.4108984555634116,-13.878683994931766,1.1520324285584622,4.4993378750758986,231.8342339859107,51.526300185224464 +53.77557228910148,-0.42197637789647496,-7.249536657710529,-0.5000497800254029,-17.797003926920357,1.185018418609399,4.2339632678902515,221.05101314891033,52.209005879982094 +54.36008937920041,-0.4967992967974098,-8.437225397893114,-0.6086852210122606,-22.376184840418,1.2252135317745365,3.910971261234011,206.94133972387158,52.91303001254382 +54.944606469299345,-0.5815938882013996,-9.728742691201678,-0.7410630528575037,-27.6806898991169,1.2741933295572765,3.5177695113374905,188.7038215862265,53.64303175010447 +55.52912355939827,-0.6765029438074703,-11.104706308706483,-0.9023722002670144,-33.76591834028289,1.3338777140987954,3.0390136043433786,165.33658856994361,54.40468852579122 +56.1136406494972,-0.7812674544568233,-12.538652143710078,-1.0989356329451032,-40.67380062540061,1.4066061841896882,2.4560055884698335,135.58358854783296,55.20491858176336 +56.69815773959613,-0.8951523346832339,-13.99937902743073,-1.3384582123317916,-48.42925380576363,1.495229538562763,1.7459610547399986,97.86487468711262,56.05215214934233 +57.282674829695054,-1.016907883547717,-15.45519893702795,-1.6303287055074418,-57.03856246502108,1.6032216210377535,0.8811160374274739,50.185428234765006,56.95666189584689 +57.86719191979398,-1.1447833032697259,-16.879887399800843,-1.9859878045565882,-66.49076342414816,1.7348154876859376,-0.17236127335316223,-9.985054995831723,57.930965590934655 diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py index 9d3bea858b..43c07deda1 100644 --- a/pvlib/pvsystem.py +++ b/pvlib/pvsystem.py @@ -20,6 +20,7 @@ from pvlib.tools import _build_kwargs from pvlib.location import Location from pvlib import irradiance, atmosphere +from pvlib import singlediode_methods # not sure if this belongs in the pvsystem module. @@ -1851,8 +1852,9 @@ def sapm_effective_irradiance(poa_direct, poa_diffuse, airmass_absolute, aoi, def singlediode(photocurrent, saturation_current, resistance_series, - resistance_shunt, nNsVth, ivcurve_pnts=None): - r''' + resistance_shunt, nNsVth, ivcurve_pnts=None, + method='lambertw'): + """ Solve the single-diode model to obtain a photovoltaic IV curve. Singlediode solves the single diode equation [1] @@ -1907,6 +1909,10 @@ def singlediode(photocurrent, saturation_current, resistance_series, Number of points in the desired IV curve. If None or 0, no IV curves will be produced. + method : str, default 'lambertw' + Determines the method used to calculate points on the IV curve. The + options are ``'lambertw'``, ``'newton'``, or ``'brentq'``. + Returns ------- OrderedDict or DataFrame @@ -1935,9 +1941,64 @@ def singlediode(photocurrent, saturation_current, resistance_series, Notes ----- - The solution employed to solve the implicit diode equation utilizes - the Lambert W function to obtain an explicit function of V=f(i) and - I=f(V) as shown in [2]. + If the method is ``'lambertw'`` then the solution employed to solve the + implicit diode equation utilizes the Lambert W function to obtain an + explicit function of :math:`V=f(I)` and :math:`I=f(V)` as shown in [2]. + + If the method is ``'newton'`` then the root-finding Newton-Raphson method + is used. It should be safe for well behaved IV-curves, but the ``'brentq'`` + method is recommended for reliability. + + If the method is ``'brentq'`` then Brent's bisection search method is used + that guarantees convergence by bounding the voltage between zero and + open-circuit. + + If the method is either ``'newton'`` or ``'brentq'`` and ``ivcurve_pnts`` + are indicated, then :func:`pvlib.singlediode_methods.bishop88` is used to + calculate the points on the IV curve points at diode voltages from zero to + open-circuit voltage with a log spacing that gets closer as voltage + increases. If the method is ``'lambertw'`` then the calculated points on + the IV curve are linearly spaced. + + The ``bishop88`` method uses an explicit solution from [4] that finds + points on the IV curve by first solving for pairs :math:`(V_d, I)` where + :math:`V_d` is the diode voltage :math:`V_d = V + I*Rs`. Then the voltage + is backed out from :math:`V_d`. Points with specific voltage, such as open + circuit, are located using the bisection search method, ``brentq``, bounded + by a zero diode voltage and an estimate of open circuit voltage given by + + .. math:: + + V_{oc, est} = n Ns V_{th} \\log \\left( \\frac{I_L}{I_0} + 1 \\right) + + We know that :math:`V_d = 0` corresponds to a voltage less than zero, and + we can also show that when :math:`V_d = V_{oc, est}`, the resulting + current is also negative, meaning that the corresponding voltage must be + in the 4th quadrant and therefore greater than the open circuit voltage + (see proof below). Therefore the entire forward-bias 1st quadrant IV-curve + is bounded, and a bisection search within these points will always find + desired condition. + + .. math:: + + I = I_L - I_0 \\left(\\exp \\left(\\frac{V_{oc, est}}{n Ns V_{th}} \\right) - 1 \\right) + - \\frac{V_{oc, est}}{R_{sh}} \\newline + + I = I_L - I_0 \\left(\\exp \\left(\\frac{n Ns V_{th} \\log \\left(\\frac{I_L}{I_0} + 1 \\right)}{n Ns V_{th}} \\right) - 1 \\right) + - \\frac{n Ns V_{th} \\log \\left(\\frac{I_L}{I_0} + 1 \\right)}{R_{sh}} \\newline + + I = I_L - I_0 \\left(\\exp \\left(\\log \\left(\\frac{I_L}{I_0} + 1 \\right) \\right) - 1 \\right) + - \\frac{n Ns V_{th} \\log \\left(\\frac{I_L}{I_0} + 1 \\right)}{R_{sh}} \\newline + + I = I_L - I_0 \\left(\\frac{I_L}{I_0} + 1 - 1 \\right) + - \\frac{n Ns V_{th} \\log \\left(\\frac{I_L}{I_0} + 1 \\right)}{R_{sh}} \\newline + + I = I_L - I_0 \\left(\\frac{I_L}{I_0} \\right) + - \\frac{n Ns V_{th} \\log \\left(\\frac{I_L}{I_0} + 1 \\right)}{R_{sh}} \\newline + + I = I_L - I_L - \\frac{n Ns V_{th} \log \\left( \\frac{I_L}{I_0} + 1 \\right)}{R_{sh}} \\newline + + I = - \\frac{n Ns V_{th} \\log \\left( \\frac{I_L}{I_0} + 1 \\right)}{R_{sh}} References ----------- @@ -1951,39 +2012,55 @@ def singlediode(photocurrent, saturation_current, resistance_series, [3] D. King et al, "Sandia Photovoltaic Array Performance Model", SAND2004-3535, Sandia National Laboratories, Albuquerque, NM + [4] "Computer simulation of the effects of electrical mismatches in + photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988) + https://doi.org/10.1016/0379-6787(88)90059-2 + See also -------- sapm calcparams_desoto - ''' - - # Compute short circuit current - i_sc = i_from_v(resistance_shunt, resistance_series, nNsVth, 0., - saturation_current, photocurrent) - - # Compute open circuit voltage - v_oc = v_from_i(resistance_shunt, resistance_series, nNsVth, 0., - saturation_current, photocurrent) - - params = {'r_sh': resistance_shunt, - 'r_s': resistance_series, - 'nNsVth': nNsVth, - 'i_0': saturation_current, - 'i_l': photocurrent} - - p_mp, v_mp = _golden_sect_DataFrame(params, 0., v_oc * 1.14, _pwr_optfcn) - - # Invert the Power-Current curve. Find the current where the inverted power - # is minimized. This is i_mp. Start the optimization at v_oc/2 - i_mp = i_from_v(resistance_shunt, resistance_series, nNsVth, v_mp, - saturation_current, photocurrent) - - # Find Ix and Ixx using Lambert W - i_x = i_from_v(resistance_shunt, resistance_series, nNsVth, 0.5 * v_oc, - saturation_current, photocurrent) - - i_xx = i_from_v(resistance_shunt, resistance_series, nNsVth, - 0.5 * (v_oc + v_mp), saturation_current, photocurrent) + pvlib.singlediode_methods.bishop88 + """ + # Calculate points on the IV curve using the LambertW solution to the + # single diode equation + if method.lower() == 'lambertw': + out = singlediode_methods._lambertw( + photocurrent, saturation_current, resistance_series, + resistance_shunt, nNsVth, ivcurve_pnts + ) + i_sc, v_oc, i_mp, v_mp, p_mp, i_x, i_xx = out[:7] + if ivcurve_pnts: + ivcurve_i, ivcurve_v = out[7:] + else: + # Calculate points on the IV curve using either 'newton' or 'brentq' + # methods. Voltages are determined by first solving the single diode + # equation for the diode voltage V_d then backing out voltage + args = (photocurrent, saturation_current, resistance_series, + resistance_shunt, nNsVth) # collect args + v_oc = singlediode_methods.bishop88_v_from_i( + 0.0, *args, method=method.lower() + ) + i_mp, v_mp, p_mp = singlediode_methods.bishop88_mpp( + *args, method=method.lower() + ) + i_sc = singlediode_methods.bishop88_i_from_v( + 0.0, *args, method=method.lower() + ) + i_x = singlediode_methods.bishop88_i_from_v( + v_oc / 2.0, *args, method=method.lower() + ) + i_xx = singlediode_methods.bishop88_i_from_v( + (v_oc + v_mp) / 2.0, *args, method=method.lower() + ) + + # calculate the IV curve if requested using bishop88 + if ivcurve_pnts: + vd = v_oc * ( + (11.0 - np.logspace(np.log10(11.0), 0.0, + ivcurve_pnts)) / 10.0 + ) + ivcurve_i, ivcurve_v, _ = singlediode_methods.bishop88(vd, *args) out = OrderedDict() out['i_sc'] = i_sc @@ -1994,13 +2071,7 @@ def singlediode(photocurrent, saturation_current, resistance_series, out['i_x'] = i_x out['i_xx'] = i_xx - # create ivcurve if ivcurve_pnts: - ivcurve_v = (np.asarray(v_oc)[..., np.newaxis] * - np.linspace(0, 1, ivcurve_pnts)) - - ivcurve_i = i_from_v(resistance_shunt, resistance_series, nNsVth, - ivcurve_v.T, saturation_current, photocurrent).T out['v'] = ivcurve_v out['i'] = ivcurve_i @@ -2011,91 +2082,57 @@ def singlediode(photocurrent, saturation_current, resistance_series, return out -# Created April,2014 -# Author: Rob Andrews, Calama Consulting - -def _golden_sect_DataFrame(params, VL, VH, func): - ''' - Vectorized golden section search for finding MPPT - from a dataframe timeseries. +def max_power_point(photocurrent, saturation_current, resistance_series, + resistance_shunt, nNsVth, method='brentq'): + """ + Given the single diode equation coefficients, calculates the maximum power + point (MPP). Parameters ---------- - params : dict - Dictionary containing scalars or arrays - of inputs to the function to be optimized. - Each row should represent an independent optimization. - - VL: float - Lower bound of the optimization - - VH: float - Upper bound of the optimization - - func: function - Function to be optimized must be in the form f(array-like, x) + photocurrent : numeric + photo-generated current [A] + saturation_current : numeric + diode reverse saturation current [A] + resistance_series : numeric + series resitance [ohms] + resistance_shunt : numeric + shunt resitance [ohms] + nNsVth : numeric + product of thermal voltage ``Vth`` [V], diode ideality factor ``n``, + and number of serices cells ``Ns`` + method : str + either ``'newton'`` or ``'brentq'`` Returns ------- - func(df,'V1') : DataFrame - function evaluated at the optimal point - - df['V1']: Dataframe - Dataframe of optimal points + OrderedDict or pandas.Datafrane + ``(i_mp, v_mp, p_mp)`` Notes ----- - This funtion will find the MAXIMUM of a function - ''' - - df = params - df['VH'] = VH - df['VL'] = VL - - err = df['VH'] - df['VL'] - errflag = True - iterations = 0 - - while errflag: - - phi = (np.sqrt(5)-1)/2*(df['VH']-df['VL']) - df['V1'] = df['VL'] + phi - df['V2'] = df['VH'] - phi - - df['f1'] = func(df, 'V1') - df['f2'] = func(df, 'V2') - df['SW_Flag'] = df['f1'] > df['f2'] - - df['VL'] = df['V2']*df['SW_Flag'] + df['VL']*(~df['SW_Flag']) - df['VH'] = df['V1']*~df['SW_Flag'] + df['VH']*(df['SW_Flag']) - - err = df['V1'] - df['V2'] - try: - errflag = (abs(err) > .01).any() - except ValueError: - errflag = (abs(err) > .01) - - iterations += 1 - - if iterations > 50: - raise Exception("EXCEPTION:iterations exeeded maximum (50)") - - return func(df, 'V1'), df['V1'] - - -def _pwr_optfcn(df, loc): - ''' - Function to find power from ``i_from_v``. - ''' - - I = i_from_v(df['r_sh'], df['r_s'], df['nNsVth'], df[loc], df['i_0'], - df['i_l']) - - return I * df[loc] + Use this function when you only want to find the maximum power point. Use + :func:`singlediode` when you need to find additional points on the IV + curve. This function uses Brent's method by default because it is + guaranteed to converge. + """ + i_mp, v_mp, p_mp = singlediode_methods.bishop88_mpp( + photocurrent, saturation_current, resistance_series, + resistance_shunt, nNsVth, method=method.lower() + ) + if isinstance(photocurrent, pd.Series): + ivp = {'i_mp': i_mp, 'v_mp': v_mp, 'p_mp': p_mp} + out = pd.DataFrame(ivp, index=photocurrent.index) + else: + out = OrderedDict() + out['i_mp'] = i_mp + out['v_mp'] = v_mp + out['p_mp'] = p_mp + return out def v_from_i(resistance_shunt, resistance_series, nNsVth, current, - saturation_current, photocurrent): + saturation_current, photocurrent, method='lambertw'): ''' Device voltage at the given device current for the single diode model. @@ -2144,6 +2181,10 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current, IV curve conditions. Often abbreviated ``I_L``. 0 <= photocurrent + method : str + Method to use: ``'lambertw'``, ``'newton'``, or ``'brentq'``. *Note*: + ``'brentq'`` is limited to 1st quadrant only. + Returns ------- current : np.ndarray or scalar @@ -2154,88 +2195,32 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current, parameters of real solar cells using Lambert W-function", Solar Energy Materials and Solar Cells, 81 (2004) 269-277. ''' - try: - from scipy.special import lambertw - except ImportError: - raise ImportError('This function requires scipy') - - # Record if inputs were all scalar - output_is_scalar = all(map(np.isscalar, - [resistance_shunt, resistance_series, nNsVth, - current, saturation_current, photocurrent])) - - # This transforms Gsh=1/Rsh, including ideal Rsh=np.inf into Gsh=0., which - # is generally more numerically stable - conductance_shunt = 1./resistance_shunt - - # Ensure that we are working with read-only views of numpy arrays - # Turns Series into arrays so that we don't have to worry about - # multidimensional broadcasting failing - Gsh, Rs, a, I, I0, IL = \ - np.broadcast_arrays(conductance_shunt, resistance_series, nNsVth, - current, saturation_current, photocurrent) - - # Intitalize output V (I might not be float64) - V = np.full_like(I, np.nan, dtype=np.float64) - - # Determine indices where 0 < Gsh requires implicit model solution - idx_p = 0. < Gsh - - # Determine indices where 0 = Gsh allows explicit model solution - idx_z = 0. == Gsh - - # Explicit solutions where Gsh=0 - if np.any(idx_z): - V[idx_z] = a[idx_z]*np.log1p((IL[idx_z] - I[idx_z])/I0[idx_z]) - \ - I[idx_z]*Rs[idx_z] - - # Only compute using LambertW if there are cases with Gsh>0 - if np.any(idx_p): - # LambertW argument, cannot be float128, may overflow to np.inf - # overflow is explicitly handled below, so ignore warnings here - with np.errstate(over='ignore'): - argW = (I0[idx_p] / (Gsh[idx_p]*a[idx_p]) * - np.exp((-I[idx_p] + IL[idx_p] + I0[idx_p]) / - (Gsh[idx_p]*a[idx_p]))) - - # lambertw typically returns complex value with zero imaginary part - # may overflow to np.inf - lambertwterm = lambertw(argW).real - - # Record indices where lambertw input overflowed output - idx_inf = np.logical_not(np.isfinite(lambertwterm)) - - # Only re-compute LambertW if it overflowed - if np.any(idx_inf): - # Calculate using log(argW) in case argW is really big - logargW = (np.log(I0[idx_p]) - np.log(Gsh[idx_p]) - - np.log(a[idx_p]) + - (-I[idx_p] + IL[idx_p] + I0[idx_p]) / - (Gsh[idx_p] * a[idx_p]))[idx_inf] - - # Three iterations of Newton-Raphson method to solve - # w+log(w)=logargW. The initial guess is w=logargW. Where direct - # evaluation (above) results in NaN from overflow, 3 iterations - # of Newton's method gives approximately 8 digits of precision. - w = logargW - for _ in range(0, 3): - w = w * (1. - np.log(w) + logargW) / (1. + w) - lambertwterm[idx_inf] = w - - # Eqn. 3 in Jain and Kapoor, 2004 - # V = -I*(Rs + Rsh) + IL*Rsh - a*lambertwterm + I0*Rsh - # Recast in terms of Gsh=1/Rsh for better numerical stability. - V[idx_p] = (IL[idx_p] + I0[idx_p] - I[idx_p])/Gsh[idx_p] - \ - I[idx_p]*Rs[idx_p] - a[idx_p]*lambertwterm - - if output_is_scalar: - return np.asscalar(V) + if method.lower() == 'lambertw': + return singlediode_methods._lambertw_v_from_i( + resistance_shunt, resistance_series, nNsVth, current, + saturation_current, photocurrent + ) else: + # Calculate points on the IV curve using either 'newton' or 'brentq' + # methods. Voltages are determined by first solving the single diode + # equation for the diode voltage V_d then backing out voltage + args = (current, photocurrent, saturation_current, + resistance_series, resistance_shunt, nNsVth) + V = singlediode_methods.bishop88_v_from_i(*args, method=method.lower()) + # find the right size and shape for returns + size, shape = singlediode_methods._get_size_and_shape(args) + if size <= 1: + if shape is not None: + V = np.tile(V, shape) + if np.isnan(V).any() and size <= 1: + V = np.repeat(V, size) + if shape is not None: + V = V.reshape(shape) return V def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage, - saturation_current, photocurrent): + saturation_current, photocurrent, method='lambertw'): ''' Device current at the given device voltage for the single diode model. @@ -2284,6 +2269,10 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage, IV curve conditions. Often abbreviated ``I_L``. 0 <= photocurrent + method : str + Method to use: ``'lambertw'``, ``'newton'``, or ``'brentq'``. *Note*: + ``'brentq'`` is limited to 1st quadrant only. + Returns ------- current : np.ndarray or scalar @@ -2294,62 +2283,27 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage, parameters of real solar cells using Lambert W-function", Solar Energy Materials and Solar Cells, 81 (2004) 269-277. ''' - try: - from scipy.special import lambertw - except ImportError: - raise ImportError('This function requires scipy') - - # Record if inputs were all scalar - output_is_scalar = all(map(np.isscalar, - [resistance_shunt, resistance_series, nNsVth, - voltage, saturation_current, photocurrent])) - - # This transforms Gsh=1/Rsh, including ideal Rsh=np.inf into Gsh=0., which - # is generally more numerically stable - conductance_shunt = 1./resistance_shunt - - # Ensure that we are working with read-only views of numpy arrays - # Turns Series into arrays so that we don't have to worry about - # multidimensional broadcasting failing - Gsh, Rs, a, V, I0, IL = \ - np.broadcast_arrays(conductance_shunt, resistance_series, nNsVth, - voltage, saturation_current, photocurrent) - - # Intitalize output I (V might not be float64) - I = np.full_like(V, np.nan, dtype=np.float64) - - # Determine indices where 0 < Rs requires implicit model solution - idx_p = 0. < Rs - - # Determine indices where 0 = Rs allows explicit model solution - idx_z = 0. == Rs - - # Explicit solutions where Rs=0 - if np.any(idx_z): - I[idx_z] = IL[idx_z] - I0[idx_z]*np.expm1(V[idx_z]/a[idx_z]) - \ - Gsh[idx_z]*V[idx_z] - - # Only compute using LambertW if there are cases with Rs>0 - # Does NOT handle possibility of overflow, github issue 298 - if np.any(idx_p): - # LambertW argument, cannot be float128, may overflow to np.inf - argW = Rs[idx_p]*I0[idx_p]/(a[idx_p]*(Rs[idx_p]*Gsh[idx_p] + 1.)) * \ - np.exp((Rs[idx_p]*(IL[idx_p] + I0[idx_p]) + V[idx_p]) / - (a[idx_p]*(Rs[idx_p]*Gsh[idx_p] + 1.))) - - # lambertw typically returns complex value with zero imaginary part - # may overflow to np.inf - lambertwterm = lambertw(argW).real - - # Eqn. 2 in Jain and Kapoor, 2004 - # I = -V/(Rs + Rsh) - (a/Rs)*lambertwterm + Rsh*(IL + I0)/(Rs + Rsh) - # Recast in terms of Gsh=1/Rsh for better numerical stability. - I[idx_p] = (IL[idx_p] + I0[idx_p] - V[idx_p]*Gsh[idx_p]) / \ - (Rs[idx_p]*Gsh[idx_p] + 1.) - (a[idx_p]/Rs[idx_p])*lambertwterm - - if output_is_scalar: - return np.asscalar(I) + if method.lower() == 'lambertw': + return singlediode_methods._lambertw_i_from_v( + resistance_shunt, resistance_series, nNsVth, voltage, + saturation_current, photocurrent + ) else: + # Calculate points on the IV curve using either 'newton' or 'brentq' + # methods. Voltages are determined by first solving the single diode + # equation for the diode voltage V_d then backing out voltage + args = (voltage, photocurrent, saturation_current, resistance_series, + resistance_shunt, nNsVth) + I = singlediode_methods.bishop88_i_from_v(*args, method=method.lower()) + # find the right size and shape for returns + size, shape = singlediode_methods._get_size_and_shape(args) + if size <= 1: + if shape is not None: + I = np.tile(I, shape) + if np.isnan(I).any() and size <= 1: + I = np.repeat(I, size) + if shape is not None: + I = I.reshape(shape) return I diff --git a/pvlib/singlediode_methods.py b/pvlib/singlediode_methods.py new file mode 100644 index 0000000000..048870d13a --- /dev/null +++ b/pvlib/singlediode_methods.py @@ -0,0 +1,550 @@ +""" +Low-level functions for solving the single diode equation. +""" + +from functools import partial +import numpy as np +from pvlib.tools import _golden_sect_DataFrame + +# Try to import brentq from scipy to use when specified in bishop88_i_from_v, +# bishop88_v_from_i, and bishop88_mpp methods below. If not imported, raises +# ImportError when 'brentq' method is specified for those methods. +try: + from scipy.optimize import brentq +except ImportError: + brentq = NotImplemented + +# FIXME: change this to newton when scipy-1.2 is released +try: + from scipy.optimize import _array_newton +except ImportError: + from pvlib.tools import _array_newton +# rename newton and set keyword arguments +newton = partial(_array_newton, tol=1e-6, maxiter=100, fprime2=None) + + +def estimate_voc(photocurrent, saturation_current, nNsVth): + """ + Rough estimate of open circuit voltage useful for bounding searches for + ``i`` of ``v`` when using :func:`~pvlib.pvsystem.singlediode`. + + Parameters + ---------- + photocurrent : numeric + photo-generated current [A] + saturation_current : numeric + diode reverse saturation current [A] + nNsVth : numeric + product of thermal voltage ``Vth`` [V], diode ideality factor ``n``, + and number of series cells ``Ns`` + + Returns + ------- + numeric + rough estimate of open circuit voltage [V] + + Notes + ----- + Calculating the open circuit voltage, :math:`V_{oc}`, of an ideal device + with infinite shunt resistance, :math:`R_{sh} \\to \\infty`, and zero + series resistance, :math:`R_s = 0`, yields the following equation [1]. As + an estimate of :math:`V_{oc}` it is useful as an upper bound for the + bisection method. + + .. math:: + + V_{oc, est}=n Ns V_{th} \\log \\left( \\frac{I_L}{I_0} + 1 \\right) + + [1] http://www.pveducation.org/pvcdrom/open-circuit-voltage + """ + + return nNsVth * np.log(np.asarray(photocurrent) / saturation_current + 1.0) + + +def bishop88(diode_voltage, photocurrent, saturation_current, + resistance_series, resistance_shunt, nNsVth, gradients=False): + """ + Explicit calculation of points on the IV curve described by the single + diode equation [1]. + + [1] "Computer simulation of the effects of electrical mismatches in + photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988) + https://doi.org/10.1016/0379-6787(88)90059-2 + + Parameters + ---------- + diode_voltage : numeric + diode voltages [V] + photocurrent : numeric + photo-generated current [A] + saturation_current : numeric + diode reverse saturation current [A] + resistance_series : numeric + series resistance [ohms] + resistance_shunt: numeric + shunt resistance [ohms] + nNsVth : numeric + product of thermal voltage ``Vth`` [V], diode ideality factor ``n``, + and number of series cells ``Ns`` + gradients : bool + False returns only I, V, and P. True also returns gradients + + Returns + ------- + tuple + currents [A], voltages [V], power [W], and optionally + :math:`\\frac{dI}{dV_d}`, :math:`\\frac{dV}{dV_d}`, + :math:`\\frac{dI}{dV}`, :math:`\\frac{dP}{dV}`, and + :math:`\\frac{d^2 P}{dV dV_d}` + """ + # calculate temporary values to simplify calculations + v_star = diode_voltage / nNsVth # non-dimensional diode voltage + g_sh = 1.0 / resistance_shunt # conductance + i = (photocurrent - saturation_current * np.expm1(v_star) + - diode_voltage * g_sh) + v = diode_voltage - i * resistance_series + retval = (i, v, i*v) + if gradients: + g_diode = saturation_current * np.exp(v_star) / nNsVth # conductance + grad_i = -g_diode - g_sh # di/dvd + grad_v = 1.0 - grad_i * resistance_series # dv/dvd + # dp/dv = d(iv)/dv = v * di/dv + i + grad = grad_i / grad_v # di/dv + grad_p = v * grad + i # dp/dv + grad2i = -g_diode / nNsVth # d2i/dvd + grad2v = -grad2i * resistance_series # d2v/dvd + grad2p = ( + grad_v * grad + v * (grad2i/grad_v - grad_i*grad2v/grad_v**2) + + grad_i + ) # d2p/dv/dvd + retval += (grad_i, grad_v, grad, grad_p, grad2p) + return retval + + +def bishop88_i_from_v(voltage, photocurrent, saturation_current, + resistance_series, resistance_shunt, nNsVth, + method='newton'): + """ + Find current given any voltage. + + Parameters + ---------- + voltage : numeric + voltage (V) in volts [V] + photocurrent : numeric + photogenerated current (Iph or IL) in amperes [A] + saturation_current : numeric + diode dark or saturation current (Io or Isat) in amperes [A] + resistance_series : numeric + series resistance (Rs) in ohms + resistance_shunt : numeric + shunt resistance (Rsh) in ohms + nNsVth : numeric + product of diode ideality factor (n), number of series cells (Ns), and + thermal voltage (Vth = k_b * T / q_e) in volts [V] + method : str + one of two optional search methods: either ``'brentq'``, a reliable and + bounded method or ``'newton'`` which is the default. + + Returns + ------- + current : numeric + current (I) at the specified voltage (V) in amperes [A] + """ + # collect args + args = (photocurrent, saturation_current, resistance_series, + resistance_shunt, nNsVth) + + def fv(x, v, *a): + # calculate voltage residual given diode voltage "x" + return bishop88(x, *a)[1] - v + + if method.lower() == 'brentq': + if brentq is NotImplemented: + raise ImportError('This function requires scipy') + # first bound the search using voc + voc_est = estimate_voc(photocurrent, saturation_current, nNsVth) + + # brentq only works with scalar inputs, so we need a set up function + # and np.vectorize to repeatedly call the optimizer with the right + # arguments for possible array input + def vd_from_brent(voc, v, iph, isat, rs, rsh, gamma): + return brentq(fv, 0.0, voc, args=(v, iph, isat, rs, rsh, gamma)) + + vd_from_brent_vectorized = np.vectorize(vd_from_brent) + vd = vd_from_brent_vectorized(voc_est, voltage, *args) + elif method.lower() == 'newton': + # make sure all args are numpy arrays if max size > 1 + # if voltage is an array, then make a copy to use for initial guess, v0 + args, v0 = _prepare_newton_inputs((voltage,), args, voltage) + vd = newton(func=lambda x, *a: fv(x, voltage, *a), x0=v0, + fprime=lambda x, *a: bishop88(x, *a, gradients=True)[4], + args=args) + else: + raise NotImplementedError("Method '%s' isn't implemented" % method) + return bishop88(vd, *args)[0] + + +def bishop88_v_from_i(current, photocurrent, saturation_current, + resistance_series, resistance_shunt, nNsVth, + method='newton'): + """ + Find voltage given any current. + + Parameters + ---------- + current : numeric + current (I) in amperes [A] + photocurrent : numeric + photogenerated current (Iph or IL) in amperes [A] + saturation_current : numeric + diode dark or saturation current (Io or Isat) in amperes [A] + resistance_series : numeric + series resistance (Rs) in ohms + resistance_shunt : numeric + shunt resistance (Rsh) in ohms + nNsVth : numeric + product of diode ideality factor (n), number of series cells (Ns), and + thermal voltage (Vth = k_b * T / q_e) in volts [V] + method : str + one of two optional search methods: either ``'brentq'``, a reliable and + bounded method or ``'newton'`` which is the default. + + Returns + ------- + voltage : numeric + voltage (V) at the specified current (I) in volts [V] + """ + # collect args + args = (photocurrent, saturation_current, resistance_series, + resistance_shunt, nNsVth) + # first bound the search using voc + voc_est = estimate_voc(photocurrent, saturation_current, nNsVth) + + def fi(x, i, *a): + # calculate current residual given diode voltage "x" + return bishop88(x, *a)[0] - i + + if method.lower() == 'brentq': + if brentq is NotImplemented: + raise ImportError('This function requires scipy') + + # brentq only works with scalar inputs, so we need a set up function + # and np.vectorize to repeatedly call the optimizer with the right + # arguments for possible array input + def vd_from_brent(voc, i, iph, isat, rs, rsh, gamma): + return brentq(fi, 0.0, voc, args=(i, iph, isat, rs, rsh, gamma)) + + vd_from_brent_vectorized = np.vectorize(vd_from_brent) + vd = vd_from_brent_vectorized(voc_est, current, *args) + elif method.lower() == 'newton': + # make sure all args are numpy arrays if max size > 1 + # if voc_est is an array, then make a copy to use for initial guess, v0 + args, v0 = _prepare_newton_inputs((current,), args, voc_est) + vd = newton(func=lambda x, *a: fi(x, current, *a), x0=v0, + fprime=lambda x, *a: bishop88(x, *a, gradients=True)[3], + args=args) + else: + raise NotImplementedError("Method '%s' isn't implemented" % method) + return bishop88(vd, *args)[1] + + +def bishop88_mpp(photocurrent, saturation_current, resistance_series, + resistance_shunt, nNsVth, method='newton'): + """ + Find max power point. + + Parameters + ---------- + photocurrent : numeric + photogenerated current (Iph or IL) in amperes [A] + saturation_current : numeric + diode dark or saturation current (Io or Isat) in amperes [A] + resistance_series : numeric + series resistance (Rs) in ohms + resistance_shunt : numeric + shunt resistance (Rsh) in ohms + nNsVth : numeric + product of diode ideality factor (n), number of series cells (Ns), and + thermal voltage (Vth = k_b * T / q_e) in volts [V] + method : str + one of two optional search methods: either ``'brentq'``, a reliable and + bounded method or ``'newton'`` which is the default. + + Returns + ------- + OrderedDict or pandas.DataFrame + max power current ``i_mp`` [A], max power voltage ``v_mp`` [V], and + max power ``p_mp`` [W] + """ + # collect args + args = (photocurrent, saturation_current, resistance_series, + resistance_shunt, nNsVth) + # first bound the search using voc + voc_est = estimate_voc(photocurrent, saturation_current, nNsVth) + + def fmpp(x, *a): + return bishop88(x, *a, gradients=True)[6] + + if method.lower() == 'brentq': + if brentq is NotImplemented: + raise ImportError('This function requires scipy') + # break out arguments for numpy.vectorize to handle broadcasting + vec_fun = np.vectorize( + lambda voc, iph, isat, rs, rsh, gamma: + brentq(fmpp, 0.0, voc, args=(iph, isat, rs, rsh, gamma)) + ) + vd = vec_fun(voc_est, *args) + elif method.lower() == 'newton': + # make sure all args are numpy arrays if max size > 1 + # if voc_est is an array, then make a copy to use for initial guess, v0 + args, v0 = _prepare_newton_inputs((), args, voc_est) + vd = newton( + func=fmpp, x0=v0, + fprime=lambda x, *a: bishop88(x, *a, gradients=True)[7], args=args + ) + else: + raise NotImplementedError("Method '%s' isn't implemented" % method) + return bishop88(vd, *args) + + +def _get_size_and_shape(args): + # find the right size and shape for returns + size, shape = 0, None # 0 or None both mean scalar + for arg in args: + try: + this_shape = arg.shape # try to get shape + except AttributeError: + this_shape = None + try: + this_size = len(arg) # try to get the size + except TypeError: + this_size = 0 + else: + this_size = arg.size # if it has shape then it also has size + if shape is None: + shape = this_shape # set the shape if None + # update size and shape + if this_size > size: + size = this_size + if this_shape is not None: + shape = this_shape + return size, shape + + +def _prepare_newton_inputs(i_or_v_tup, args, v0): + # broadcast arguments for newton method + # the first argument should be a tuple, eg: (i,), (v,) or () + size, shape = _get_size_and_shape(i_or_v_tup + args) + if size > 1: + args = [np.asarray(arg) for arg in args] + # newton uses initial guess for the output shape + # copy v0 to a new array and broadcast it to the shape of max size + if shape is not None: + v0 = np.broadcast_to(v0, shape).copy() + return args, v0 + + +def _lambertw_v_from_i(resistance_shunt, resistance_series, nNsVth, current, + saturation_current, photocurrent): + try: + from scipy.special import lambertw + except ImportError: + raise ImportError('This function requires scipy') + + # Record if inputs were all scalar + output_is_scalar = all(map(np.isscalar, + [resistance_shunt, resistance_series, nNsVth, + current, saturation_current, photocurrent])) + + # This transforms Gsh=1/Rsh, including ideal Rsh=np.inf into Gsh=0., which + # is generally more numerically stable + conductance_shunt = 1. / resistance_shunt + + # Ensure that we are working with read-only views of numpy arrays + # Turns Series into arrays so that we don't have to worry about + # multidimensional broadcasting failing + Gsh, Rs, a, I, I0, IL = \ + np.broadcast_arrays(conductance_shunt, resistance_series, nNsVth, + current, saturation_current, photocurrent) + + # Intitalize output V (I might not be float64) + V = np.full_like(I, np.nan, dtype=np.float64) + + # Determine indices where 0 < Gsh requires implicit model solution + idx_p = 0. < Gsh + + # Determine indices where 0 = Gsh allows explicit model solution + idx_z = 0. == Gsh + + # Explicit solutions where Gsh=0 + if np.any(idx_z): + V[idx_z] = a[idx_z] * np.log1p((IL[idx_z] - I[idx_z]) / I0[idx_z]) - \ + I[idx_z] * Rs[idx_z] + + # Only compute using LambertW if there are cases with Gsh>0 + if np.any(idx_p): + # LambertW argument, cannot be float128, may overflow to np.inf + # overflow is explicitly handled below, so ignore warnings here + with np.errstate(over='ignore'): + argW = (I0[idx_p] / (Gsh[idx_p] * a[idx_p]) * + np.exp((-I[idx_p] + IL[idx_p] + I0[idx_p]) / + (Gsh[idx_p] * a[idx_p]))) + + # lambertw typically returns complex value with zero imaginary part + # may overflow to np.inf + lambertwterm = lambertw(argW).real + + # Record indices where lambertw input overflowed output + idx_inf = np.logical_not(np.isfinite(lambertwterm)) + + # Only re-compute LambertW if it overflowed + if np.any(idx_inf): + # Calculate using log(argW) in case argW is really big + logargW = (np.log(I0[idx_p]) - np.log(Gsh[idx_p]) - + np.log(a[idx_p]) + + (-I[idx_p] + IL[idx_p] + I0[idx_p]) / + (Gsh[idx_p] * a[idx_p]))[idx_inf] + + # Three iterations of Newton-Raphson method to solve + # w+log(w)=logargW. The initial guess is w=logargW. Where direct + # evaluation (above) results in NaN from overflow, 3 iterations + # of Newton's method gives approximately 8 digits of precision. + w = logargW + for _ in range(0, 3): + w = w * (1. - np.log(w) + logargW) / (1. + w) + lambertwterm[idx_inf] = w + + # Eqn. 3 in Jain and Kapoor, 2004 + # V = -I*(Rs + Rsh) + IL*Rsh - a*lambertwterm + I0*Rsh + # Recast in terms of Gsh=1/Rsh for better numerical stability. + V[idx_p] = (IL[idx_p] + I0[idx_p] - I[idx_p]) / Gsh[idx_p] - \ + I[idx_p] * Rs[idx_p] - a[idx_p] * lambertwterm + + if output_is_scalar: + return np.asscalar(V) + else: + return V + + +def _lambertw_i_from_v(resistance_shunt, resistance_series, nNsVth, voltage, + saturation_current, photocurrent): + try: + from scipy.special import lambertw + except ImportError: + raise ImportError('This function requires scipy') + + # Record if inputs were all scalar + output_is_scalar = all(map(np.isscalar, + [resistance_shunt, resistance_series, nNsVth, + voltage, saturation_current, photocurrent])) + + # This transforms Gsh=1/Rsh, including ideal Rsh=np.inf into Gsh=0., which + # is generally more numerically stable + conductance_shunt = 1. / resistance_shunt + + # Ensure that we are working with read-only views of numpy arrays + # Turns Series into arrays so that we don't have to worry about + # multidimensional broadcasting failing + Gsh, Rs, a, V, I0, IL = \ + np.broadcast_arrays(conductance_shunt, resistance_series, nNsVth, + voltage, saturation_current, photocurrent) + + # Intitalize output I (V might not be float64) + I = np.full_like(V, np.nan, dtype=np.float64) + + # Determine indices where 0 < Rs requires implicit model solution + idx_p = 0. < Rs + + # Determine indices where 0 = Rs allows explicit model solution + idx_z = 0. == Rs + + # Explicit solutions where Rs=0 + if np.any(idx_z): + I[idx_z] = IL[idx_z] - I0[idx_z] * np.expm1(V[idx_z] / a[idx_z]) - \ + Gsh[idx_z] * V[idx_z] + + # Only compute using LambertW if there are cases with Rs>0 + # Does NOT handle possibility of overflow, github issue 298 + if np.any(idx_p): + # LambertW argument, cannot be float128, may overflow to np.inf + argW = Rs[idx_p] * I0[idx_p] / ( + a[idx_p] * (Rs[idx_p] * Gsh[idx_p] + 1.)) * \ + np.exp((Rs[idx_p] * (IL[idx_p] + I0[idx_p]) + V[idx_p]) / + (a[idx_p] * (Rs[idx_p] * Gsh[idx_p] + 1.))) + + # lambertw typically returns complex value with zero imaginary part + # may overflow to np.inf + lambertwterm = lambertw(argW).real + + # Eqn. 2 in Jain and Kapoor, 2004 + # I = -V/(Rs + Rsh) - (a/Rs)*lambertwterm + Rsh*(IL + I0)/(Rs + Rsh) + # Recast in terms of Gsh=1/Rsh for better numerical stability. + I[idx_p] = (IL[idx_p] + I0[idx_p] - V[idx_p] * Gsh[idx_p]) / \ + (Rs[idx_p] * Gsh[idx_p] + 1.) - ( + a[idx_p] / Rs[idx_p]) * lambertwterm + + if output_is_scalar: + return np.asscalar(I) + else: + return I + + +def _lambertw(photocurrent, saturation_current, resistance_series, + resistance_shunt, nNsVth, ivcurve_pnts=None): + # Compute short circuit current + i_sc = _lambertw_i_from_v(resistance_shunt, resistance_series, nNsVth, 0., + saturation_current, photocurrent) + + # Compute open circuit voltage + v_oc = _lambertw_v_from_i(resistance_shunt, resistance_series, nNsVth, 0., + saturation_current, photocurrent) + + params = {'r_sh': resistance_shunt, + 'r_s': resistance_series, + 'nNsVth': nNsVth, + 'i_0': saturation_current, + 'i_l': photocurrent} + + # Find the voltage, v_mp, where the power is maximized. + # Start the golden section search at v_oc * 1.14 + p_mp, v_mp = _golden_sect_DataFrame(params, 0., v_oc * 1.14, + _pwr_optfcn) + + # Find Imp using Lambert W + i_mp = _lambertw_i_from_v(resistance_shunt, resistance_series, nNsVth, + v_mp, saturation_current, photocurrent) + + # Find Ix and Ixx using Lambert W + i_x = _lambertw_i_from_v(resistance_shunt, resistance_series, nNsVth, + 0.5 * v_oc, saturation_current, photocurrent) + + i_xx = _lambertw_i_from_v(resistance_shunt, resistance_series, nNsVth, + 0.5 * (v_oc + v_mp), saturation_current, + photocurrent) + + out = (i_sc, v_oc, i_mp, v_mp, p_mp, i_x, i_xx) + + # create ivcurve + if ivcurve_pnts: + ivcurve_v = (np.asarray(v_oc)[..., np.newaxis] * + np.linspace(0, 1, ivcurve_pnts)) + + ivcurve_i = _lambertw_i_from_v(resistance_shunt, resistance_series, + nNsVth, ivcurve_v.T, saturation_current, + photocurrent).T + + out += (ivcurve_i, ivcurve_v) + + return out + + +def _pwr_optfcn(df, loc): + ''' + Function to find power from ``i_from_v``. + ''' + + I = _lambertw_i_from_v(df['r_sh'], df['r_s'], df['nNsVth'], df[loc], + df['i_0'], df['i_l']) + + return I * df[loc] diff --git a/pvlib/test/test_numerical_precision.py b/pvlib/test/test_numerical_precision.py new file mode 100644 index 0000000000..8d4ebd3d7d --- /dev/null +++ b/pvlib/test/test_numerical_precision.py @@ -0,0 +1,122 @@ +""" +Test numerical precision of explicit single diode calculation using symbolic +mathematics. SymPy is a computer algebra system, that uses infinite precision +symbols instead of standard floating point and integer computer number types. +http://docs.sympy.org/latest/modules/evalf.html#accuracy-and-error-handling + +This module can be executed from the command line to generate a high precision +dataset of I-V curve points to test the explicit single diode calculations +:func:`pvlib.singlediode_methods.bishop88`:: + + $ python test_numeric_precision.py + +This generates a file in the pvlib data folder, which is specified by the +constant ``DATA_PATH``. When the test is run using ``pytest`` it will compare +the values calculated by :func:`pvlib.singlediode_methods.bishop88` with the +high-precision values generated with SymPy. +""" + +import logging +import os +import numpy as np +import pandas as pd +from pvlib import pvsystem +from pvlib.singlediode_methods import bishop88, estimate_voc + +logging.basicConfig() +LOGGER = logging.getLogger(__name__) +LOGGER.setLevel(logging.DEBUG) +TEST_DATA = 'bishop88_numerical_precision.csv' +TEST_PATH = os.path.dirname(os.path.abspath(__file__)) +PVLIB_PATH = os.path.dirname(TEST_PATH) +DATA_PATH = os.path.join(PVLIB_PATH, 'data', TEST_DATA) +POA = 888 +TCELL = 55 +CECMOD = pvsystem.retrieve_sam('cecmod') +# get module from cecmod and apply temp/irrad desoto corrections +SPR_E20_327 = CECMOD.SunPower_SPR_E20_327 +ARGS = pvsystem.calcparams_desoto( + effective_irradiance=POA, temp_cell=TCELL, + alpha_sc=SPR_E20_327.alpha_sc, a_ref=SPR_E20_327.a_ref, + I_L_ref=SPR_E20_327.I_L_ref, I_o_ref=SPR_E20_327.I_o_ref, + R_sh_ref=SPR_E20_327.R_sh_ref, R_s=SPR_E20_327.R_s, + EgRef=1.121, dEgdT=-0.0002677 +) +IL, I0, RS, RSH, NNSVTH = ARGS +IVCURVE_NPTS = 100 + +try: + from sympy import symbols, exp as sy_exp +except ImportError as exc: + LOGGER.exception(exc) + symbols = NotImplemented + sy_exp = NotImplemented + + +def generate_numerical_precision(): + """ + Generate expected data with infinite numerical precision using SymPy. + :return: dataframe of expected values + """ + if symbols is NotImplemented: + LOGGER.critical("SymPy is required to generate expected data.") + raise ImportError("could not import sympy") + il, io, rs, rsh, nnsvt, vd = symbols('il, io, rs, rsh, nnsvt, vd') + a = sy_exp(vd / nnsvt) + b = 1.0 / rsh + i = il - io * (a - 1.0) - vd * b + v = vd - i * rs + c = io * a / nnsvt + grad_i = - c - b # di/dvd + grad_v = 1.0 - grad_i * rs # dv/dvd + # dp/dv = d(iv)/dv = v * di/dv + i + grad = grad_i / grad_v # di/dv + p = i * v + grad_p = v * grad + i # dp/dv + grad2i = -c / nnsvt + grad2v = -grad2i * rs + grad2p = ( + grad_v * grad + v * (grad2i/grad_v - grad_i*grad2v/grad_v**2) + grad_i + ) + # generate exact values + data = dict(zip((il, io, rs, rsh, nnsvt), ARGS)) + vdtest = np.linspace(0, estimate_voc(IL, I0, NNSVTH), IVCURVE_NPTS) + expected = [] + for test in vdtest: + data[vd] = test + test_data = { + 'i': np.float64(i.evalf(subs=data)), + 'v': np.float64(v.evalf(subs=data)), + 'p': np.float64(p.evalf(subs=data)), + 'grad_i': np.float64(grad_i.evalf(subs=data)), + 'grad_v': np.float64(grad_v.evalf(subs=data)), + 'grad': np.float64(grad.evalf(subs=data)), + 'grad_p': np.float64(grad_p.evalf(subs=data)), + 'grad2p': np.float64(grad2p.evalf(subs=data)) + } + LOGGER.debug(test_data) + expected.append(test_data) + return pd.DataFrame(expected, index=vdtest) + + +def test_numerical_precision(): + """ + Test that there are no numerical errors due to floating point arithmetic. + """ + expected = pd.read_csv(DATA_PATH) + vdtest = np.linspace(0, estimate_voc(IL, I0, NNSVTH), IVCURVE_NPTS) + results = bishop88(vdtest, *ARGS, gradients=True) + assert np.allclose(expected['i'], results[0]) + assert np.allclose(expected['v'], results[1]) + assert np.allclose(expected['p'], results[2]) + assert np.allclose(expected['grad_i'], results[3]) + assert np.allclose(expected['grad_v'], results[4]) + assert np.allclose(expected['grad'], results[5]) + assert np.allclose(expected['grad_p'], results[6]) + assert np.allclose(expected['grad2p'], results[7]) + + +if __name__ == '__main__': + expected = generate_numerical_precision() + expected.to_csv(DATA_PATH) + test_numerical_precision() diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py index 6931543824..729a654ccc 100644 --- a/pvlib/test/test_pvsystem.py +++ b/pvlib/test/test_pvsystem.py @@ -608,7 +608,10 @@ def fixture_v_from_i(request): @requires_scipy -def test_v_from_i(fixture_v_from_i): +@pytest.mark.parametrize( + 'method, atol', [('lambertw', 1e-11), ('brentq', 1e-11), ('newton', 1e-8)] +) +def test_v_from_i(fixture_v_from_i, method, atol): # Solution set loaded from fixture Rsh = fixture_v_from_i['Rsh'] Rs = fixture_v_from_i['Rs'] @@ -618,10 +621,7 @@ def test_v_from_i(fixture_v_from_i): IL = fixture_v_from_i['IL'] V_expected = fixture_v_from_i['V_expected'] - # Convergence criteria - atol = 1.e-11 - - V = pvsystem.v_from_i(Rsh, Rs, nNsVth, I, I0, IL) + V = pvsystem.v_from_i(Rsh, Rs, nNsVth, I, I0, IL, method=method) assert(isinstance(V, type(V_expected))) if isinstance(V, type(np.ndarray)): assert(isinstance(V.dtype, type(V_expected.dtype))) @@ -629,43 +629,68 @@ def test_v_from_i(fixture_v_from_i): assert_allclose(V, V_expected, atol=atol) +@requires_scipy +def test_i_from_v_from_i(fixture_v_from_i): + # Solution set loaded from fixture + Rsh = fixture_v_from_i['Rsh'] + Rs = fixture_v_from_i['Rs'] + nNsVth = fixture_v_from_i['nNsVth'] + I = fixture_v_from_i['I'] + I0 = fixture_v_from_i['I0'] + IL = fixture_v_from_i['IL'] + V = fixture_v_from_i['V_expected'] + + # Convergence criteria + atol = 1.e-11 + + I_expected = pvsystem.i_from_v(Rsh, Rs, nNsVth, V, I0, IL, + method='lambertw') + assert_allclose(I, I_expected, atol=atol) + I = pvsystem.i_from_v(Rsh, Rs, nNsVth, V, I0, IL) + assert(isinstance(I, type(I_expected))) + if isinstance(I, type(np.ndarray)): + assert(isinstance(I.dtype, type(I_expected.dtype))) + assert(I.shape == I_expected.shape) + assert_allclose(I, I_expected, atol=atol) + + @pytest.fixture(params=[ { # Can handle all python scalar inputs 'Rsh': 20., 'Rs': 0.1, 'nNsVth': 0.5, - 'V': 40., + 'V': 7.5049875193450521, 'I0': 6.e-7, 'IL': 7., - 'I_expected': -299.746389916 + 'I_expected': 3. }, { # Can handle all rank-0 array inputs 'Rsh': np.array(20.), 'Rs': np.array(0.1), 'nNsVth': np.array(0.5), - 'V': np.array(40.), + 'V': np.array(7.5049875193450521), 'I0': np.array(6.e-7), 'IL': np.array(7.), - 'I_expected': np.array(-299.746389916) + 'I_expected': np.array(3.) }, { # Can handle all rank-1 singleton array inputs 'Rsh': np.array([20.]), 'Rs': np.array([0.1]), 'nNsVth': np.array([0.5]), - 'V': np.array([40.]), + 'V': np.array([7.5049875193450521]), 'I0': np.array([6.e-7]), 'IL': np.array([7.]), - 'I_expected': np.array([-299.746389916]) + 'I_expected': np.array([3.]) }, { # Can handle all rank-1 non-singleton array inputs with a zero # series resistance, Rs=0 gives I=IL=Isc at V=0 'Rsh': np.array([20., 20.]), 'Rs': np.array([0., 0.1]), 'nNsVth': np.array([0.5, 0.5]), - 'V': np.array([0., 40.]), + 'V': np.array([0., 7.5049875193450521]), 'I0': np.array([6.e-7, 6.e-7]), 'IL': np.array([7., 7.]), - 'I_expected': np.array([7., -299.746389916]) + 'I_expected': np.array([7., 3.]) }, { # Can handle mixed inputs with a rank-2 array with zero series # resistance, Rs=0 gives I=IL=Isc at V=0 @@ -693,17 +718,20 @@ def test_v_from_i(fixture_v_from_i): 'Rsh': np.inf, 'Rs': 0.1, 'nNsVth': 0.5, - 'V': 40., + 'V': 7.5049875193450521, 'I0': 6.e-7, 'IL': 7., - 'I_expected': -299.7383436645412 + 'I_expected': 3.2244873645510923 }]) def fixture_i_from_v(request): return request.param @requires_scipy -def test_i_from_v(fixture_i_from_v): +@pytest.mark.parametrize( + 'method, atol', [('lambertw', 1e-11), ('brentq', 1e-11), ('newton', 1e-11)] +) +def test_i_from_v(fixture_i_from_v, method, atol): # Solution set loaded from fixture Rsh = fixture_i_from_v['Rsh'] Rs = fixture_i_from_v['Rs'] @@ -713,10 +741,7 @@ def test_i_from_v(fixture_i_from_v): IL = fixture_i_from_v['IL'] I_expected = fixture_i_from_v['I_expected'] - # Convergence criteria - atol = 1.e-11 - - I = pvsystem.i_from_v(Rsh, Rs, nNsVth, V, I0, IL) + I = pvsystem.i_from_v(Rsh, Rs, nNsVth, V, I0, IL, method=method) assert(isinstance(I, type(I_expected))) if isinstance(I, type(np.ndarray)): assert(isinstance(I.dtype, type(I_expected.dtype))) @@ -728,26 +753,102 @@ def test_i_from_v(fixture_i_from_v): def test_PVSystem_i_from_v(mocker): system = pvsystem.PVSystem() m = mocker.patch('pvlib.pvsystem.i_from_v', autospec=True) - args = (20, .1, .5, 40, 6e-7, 7) + args = (20, 0.1, 0.5, 7.5049875193450521, 6e-7, 7) system.i_from_v(*args) m.assert_called_once_with(*args) +@requires_scipy +def test_i_from_v_size(): + with pytest.raises(ValueError): + pvsystem.i_from_v(20, [0.1] * 2, 0.5, [7.5] * 3, 6.0e-7, 7.0) + with pytest.raises(ValueError): + pvsystem.i_from_v(20, [0.1] * 2, 0.5, [7.5] * 3, 6.0e-7, 7.0, + method='brentq') + with pytest.raises(ValueError): + pvsystem.i_from_v(20, 0.1, 0.5, [7.5] * 3, 6.0e-7, np.array([7., 7.]), + method='newton') + + +@requires_scipy +def test_v_from_i_size(): + with pytest.raises(ValueError): + pvsystem.v_from_i(20, [0.1] * 2, 0.5, [3.0] * 3, 6.0e-7, 7.0) + with pytest.raises(ValueError): + pvsystem.v_from_i(20, [0.1] * 2, 0.5, [3.0] * 3, 6.0e-7, 7.0, + method='brentq') + with pytest.raises(ValueError): + pvsystem.v_from_i(20, [0.1], 0.5, [3.0] * 3, 6.0e-7, np.array([7., 7.]), + method='newton') + + +@requires_scipy +def test_mpp_floats(): + """test max_power_point""" + IL, I0, Rs, Rsh, nNsVth = (7, 6e-7, .1, 20, .5) + out = pvsystem.max_power_point(IL, I0, Rs, Rsh, nNsVth, method='brentq') + expected = {'i_mp': 6.1362673597376753, # 6.1390251797935704, lambertw + 'v_mp': 6.2243393757884284, # 6.221535886625464, lambertw + 'p_mp': 38.194210547580511} # 38.194165464983037} lambertw + assert isinstance(out, dict) + for k, v in out.items(): + assert np.isclose(v, expected[k]) + out = pvsystem.max_power_point(IL, I0, Rs, Rsh, nNsVth, method='newton') + for k, v in out.items(): + assert np.isclose(v, expected[k]) + + +@requires_scipy +def test_mpp_array(): + """test max_power_point""" + IL, I0, Rs, Rsh, nNsVth = (np.array([7, 7]), 6e-7, .1, 20, .5) + out = pvsystem.max_power_point(IL, I0, Rs, Rsh, nNsVth, method='brentq') + expected = {'i_mp': [6.1362673597376753] * 2, + 'v_mp': [6.2243393757884284] * 2, + 'p_mp': [38.194210547580511] * 2} + assert isinstance(out, dict) + for k, v in out.items(): + assert np.allclose(v, expected[k]) + out = pvsystem.max_power_point(IL, I0, Rs, Rsh, nNsVth, method='newton') + for k, v in out.items(): + assert np.allclose(v, expected[k]) + + +@requires_scipy +def test_mpp_series(): + """test max_power_point""" + idx = ['2008-02-17T11:30:00-0800', '2008-02-17T12:30:00-0800'] + IL, I0, Rs, Rsh, nNsVth = (np.array([7, 7]), 6e-7, .1, 20, .5) + IL = pd.Series(IL, index=idx) + out = pvsystem.max_power_point(IL, I0, Rs, Rsh, nNsVth, method='brentq') + expected = pd.DataFrame({'i_mp': [6.1362673597376753] * 2, + 'v_mp': [6.2243393757884284] * 2, + 'p_mp': [38.194210547580511] * 2}, + index=idx) + assert isinstance(out, pd.DataFrame) + for k, v in out.items(): + assert np.allclose(v, expected[k]) + out = pvsystem.max_power_point(IL, I0, Rs, Rsh, nNsVth, method='newton') + for k, v in out.items(): + assert np.allclose(v, expected[k]) + + @requires_scipy def test_singlediode_series(cec_module_params): times = pd.DatetimeIndex(start='2015-01-01', periods=2, freq='12H') effective_irradiance = pd.Series([0.0, 800.0], index=times) IL, I0, Rs, Rsh, nNsVth = pvsystem.calcparams_desoto( - effective_irradiance, - temp_cell=25, - alpha_sc=cec_module_params['alpha_sc'], - a_ref=cec_module_params['a_ref'], - I_L_ref=cec_module_params['I_L_ref'], - I_o_ref=cec_module_params['I_o_ref'], - R_sh_ref=cec_module_params['R_sh_ref'], - R_s=cec_module_params['R_s'], - EgRef=1.121, - dEgdT=-0.0002677) + effective_irradiance, + temp_cell=25, + alpha_sc=cec_module_params['alpha_sc'], + a_ref=cec_module_params['a_ref'], + I_L_ref=cec_module_params['I_L_ref'], + I_o_ref=cec_module_params['I_o_ref'], + R_sh_ref=cec_module_params['R_sh_ref'], + R_s=cec_module_params['R_s'], + EgRef=1.121, + dEgdT=-0.0002677 + ) out = pvsystem.singlediode(IL, I0, Rs, Rsh, nNsVth) assert isinstance(out, pd.DataFrame) @@ -762,7 +863,8 @@ def test_singlediode_array(): saturation_current = 1.943e-09 sd = pvsystem.singlediode(photocurrent, saturation_current, - resistance_series, resistance_shunt, nNsVth) + resistance_series, resistance_shunt, nNsVth, + method='lambertw') expected = np.array([ 0. , 0.54538398, 1.43273966, 2.36328163, 3.29255606, @@ -771,12 +873,21 @@ def test_singlediode_array(): assert_allclose(sd['i_mp'], expected, atol=0.01) + sd = pvsystem.singlediode(photocurrent, saturation_current, + resistance_series, resistance_shunt, nNsVth) + + expected = pvsystem.i_from_v(resistance_shunt, resistance_series, nNsVth, + sd['v_mp'], saturation_current, photocurrent, + method='lambertw') + + assert_allclose(sd['i_mp'], expected, atol=0.01) + @requires_scipy def test_singlediode_floats(sam_data): module = 'Example_Module' module_parameters = sam_data['cecmod'][module] - out = pvsystem.singlediode(7, 6e-7, .1, 20, .5) + out = pvsystem.singlediode(7, 6e-7, .1, 20, .5, method='lambertw') expected = {'i_xx': 4.2498, 'i_mp': 6.1275, 'v_oc': 8.1063, @@ -796,7 +907,7 @@ def test_singlediode_floats(sam_data): @requires_scipy def test_singlediode_floats_ivcurve(): - out = pvsystem.singlediode(7, 6e-7, .1, 20, .5, ivcurve_pnts=3) + out = pvsystem.singlediode(7, 6e-7, .1, 20, .5, ivcurve_pnts=3, method='lambertw') expected = {'i_xx': 4.2498, 'i_mp': 6.1275, 'v_oc': 8.1063, @@ -827,7 +938,8 @@ def test_singlediode_series_ivcurve(cec_module_params): EgRef=1.121, dEgdT=-0.0002677) - out = pvsystem.singlediode(IL, I0, Rs, Rsh, nNsVth, ivcurve_pnts=3) + out = pvsystem.singlediode(IL, I0, Rs, Rsh, nNsVth, ivcurve_pnts=3, + method='lambertw') expected = OrderedDict([('i_sc', array([0., 3.01054475, 6.00675648])), ('v_oc', array([0., 9.96886962, 10.29530483])), @@ -848,6 +960,20 @@ def test_singlediode_series_ivcurve(cec_module_params): for k, v in out.items(): assert_allclose(v, expected[k], atol=1e-2) + out = pvsystem.singlediode(IL, I0, Rs, Rsh, nNsVth, ivcurve_pnts=3) + + expected['i_mp'] = pvsystem.i_from_v(Rsh, Rs, nNsVth, out['v_mp'], I0, IL, + method='lambertw') + expected['v_mp'] = pvsystem.v_from_i(Rsh, Rs, nNsVth, out['i_mp'], I0, IL, + method='lambertw') + expected['i'] = pvsystem.i_from_v(Rsh, Rs, nNsVth, out['v'].T, I0, IL, + method='lambertw').T + expected['v'] = pvsystem.v_from_i(Rsh, Rs, nNsVth, out['i'].T, I0, IL, + method='lambertw').T + + for k, v in out.items(): + assert_allclose(v, expected[k], atol=1e-2) + def test_scale_voltage_current_power(sam_data): data = pd.DataFrame( diff --git a/pvlib/test/test_singlediode_methods.py b/pvlib/test/test_singlediode_methods.py new file mode 100644 index 0000000000..93ad8020ed --- /dev/null +++ b/pvlib/test/test_singlediode_methods.py @@ -0,0 +1,127 @@ +""" +testing single-diode methods using JW Bishop 1988 +""" + +import numpy as np +from pvlib import pvsystem +from conftest import requires_scipy + +POA = 888 +TCELL = 55 +CECMOD = pvsystem.retrieve_sam('cecmod') + + +@requires_scipy +def test_newton_spr_e20_327(): + """test pvsystem.singlediode with Newton method on SPR-E20-327""" + spr_e20_327 = CECMOD.SunPower_SPR_E20_327 + x = pvsystem.calcparams_desoto( + effective_irradiance=POA, temp_cell=TCELL, + alpha_sc=spr_e20_327.alpha_sc, a_ref=spr_e20_327.a_ref, + I_L_ref=spr_e20_327.I_L_ref, I_o_ref=spr_e20_327.I_o_ref, + R_sh_ref=spr_e20_327.R_sh_ref, R_s=spr_e20_327.R_s, + EgRef=1.121, dEgdT=-0.0002677) + il, io, rs, rsh, nnsvt = x + pvs = pvsystem.singlediode(*x, method='lambertw') + out = pvsystem.singlediode(*x, method='newton') + isc, voc, imp, vmp, pmp, ix, ixx = out.values() + assert np.isclose(pvs['i_sc'], isc) + assert np.isclose(pvs['v_oc'], voc) + # the singlediode method doesn't actually get the MPP correct + pvs_imp = pvsystem.i_from_v(rsh, rs, nnsvt, vmp, io, il, method='lambertw') + pvs_vmp = pvsystem.v_from_i(rsh, rs, nnsvt, imp, io, il, method='lambertw') + assert np.isclose(pvs_imp, imp) + assert np.isclose(pvs_vmp, vmp) + assert np.isclose(pvs['p_mp'], pmp) + assert np.isclose(pvs['i_x'], ix) + pvs_ixx = pvsystem.i_from_v(rsh, rs, nnsvt, (voc + vmp)/2, io, il, + method='lambertw') + assert np.isclose(pvs_ixx, ixx) + return isc, voc, imp, vmp, pmp, pvs + + +@requires_scipy +def test_newton_fs_495(): + """test pvsystem.singlediode with Newton method on FS495""" + fs_495 = CECMOD.First_Solar_FS_495 + x = pvsystem.calcparams_desoto( + effective_irradiance=POA, temp_cell=TCELL, + alpha_sc=fs_495.alpha_sc, a_ref=fs_495.a_ref, I_L_ref=fs_495.I_L_ref, + I_o_ref=fs_495.I_o_ref, R_sh_ref=fs_495.R_sh_ref, R_s=fs_495.R_s, + EgRef=1.475, dEgdT=-0.0003) + il, io, rs, rsh, nnsvt = x + x += (101, ) + pvs = pvsystem.singlediode(*x, method='lambertw') + out = pvsystem.singlediode(*x, method='newton') + isc, voc, imp, vmp, pmp, ix, ixx, i, v = out.values() + assert np.isclose(pvs['i_sc'], isc) + assert np.isclose(pvs['v_oc'], voc) + # the singlediode method doesn't actually get the MPP correct + pvs_imp = pvsystem.i_from_v(rsh, rs, nnsvt, vmp, io, il, method='lambertw') + pvs_vmp = pvsystem.v_from_i(rsh, rs, nnsvt, imp, io, il, method='lambertw') + assert np.isclose(pvs_imp, imp) + assert np.isclose(pvs_vmp, vmp) + assert np.isclose(pvs['p_mp'], pmp) + assert np.isclose(pvs['i_x'], ix) + pvs_ixx = pvsystem.i_from_v(rsh, rs, nnsvt, (voc + vmp)/2, io, il, + method='lambertw') + assert np.isclose(pvs_ixx, ixx) + return isc, voc, imp, vmp, pmp, i, v, pvs + + +@requires_scipy +def test_brentq_spr_e20_327(): + """test pvsystem.singlediode with Brent method on SPR-E20-327""" + spr_e20_327 = CECMOD.SunPower_SPR_E20_327 + x = pvsystem.calcparams_desoto( + effective_irradiance=POA, temp_cell=TCELL, + alpha_sc=spr_e20_327.alpha_sc, a_ref=spr_e20_327.a_ref, + I_L_ref=spr_e20_327.I_L_ref, I_o_ref=spr_e20_327.I_o_ref, + R_sh_ref=spr_e20_327.R_sh_ref, R_s=spr_e20_327.R_s, + EgRef=1.121, dEgdT=-0.0002677) + il, io, rs, rsh, nnsvt = x + pvs = pvsystem.singlediode(*x, method='lambertw') + out = pvsystem.singlediode(*x, method='brentq') + isc, voc, imp, vmp, pmp, ix, ixx = out.values() + assert np.isclose(pvs['i_sc'], isc) + assert np.isclose(pvs['v_oc'], voc) + # the singlediode method doesn't actually get the MPP correct + pvs_imp = pvsystem.i_from_v(rsh, rs, nnsvt, vmp, io, il, method='lambertw') + pvs_vmp = pvsystem.v_from_i(rsh, rs, nnsvt, imp, io, il, method='lambertw') + assert np.isclose(pvs_imp, imp) + assert np.isclose(pvs_vmp, vmp) + assert np.isclose(pvs['p_mp'], pmp) + assert np.isclose(pvs['i_x'], ix) + pvs_ixx = pvsystem.i_from_v(rsh, rs, nnsvt, (voc + vmp)/2, io, il, + method='lambertw') + assert np.isclose(pvs_ixx, ixx) + return isc, voc, imp, vmp, pmp, pvs + + +@requires_scipy +def test_brentq_fs_495(): + """test pvsystem.singlediode with Brent method on SPR-E20-327""" + fs_495 = CECMOD.First_Solar_FS_495 + x = pvsystem.calcparams_desoto( + effective_irradiance=POA, temp_cell=TCELL, + alpha_sc=fs_495.alpha_sc, a_ref=fs_495.a_ref, I_L_ref=fs_495.I_L_ref, + I_o_ref=fs_495.I_o_ref, R_sh_ref=fs_495.R_sh_ref, R_s=fs_495.R_s, + EgRef=1.475, dEgdT=-0.0003) + il, io, rs, rsh, nnsvt = x + x += (101, ) + pvs = pvsystem.singlediode(*x, method='lambertw') + out = pvsystem.singlediode(*x, method='brentq') + isc, voc, imp, vmp, pmp, ix, ixx, i, v = out.values() + assert np.isclose(pvs['i_sc'], isc) + assert np.isclose(pvs['v_oc'], voc) + # the singlediode method doesn't actually get the MPP correct + pvs_imp = pvsystem.i_from_v(rsh, rs, nnsvt, vmp, io, il, method='lambertw') + pvs_vmp = pvsystem.v_from_i(rsh, rs, nnsvt, imp, io, il, method='lambertw') + assert np.isclose(pvs_imp, imp) + assert np.isclose(pvs_vmp, vmp) + assert np.isclose(pvs['p_mp'], pmp) + assert np.isclose(pvs['i_x'], ix) + pvs_ixx = pvsystem.i_from_v(rsh, rs, nnsvt, (voc + vmp)/2, io, il, + method='lambertw') + assert np.isclose(pvs_ixx, ixx) + return isc, voc, imp, vmp, pmp, i, v, pvs diff --git a/pvlib/tools.py b/pvlib/tools.py index a722500f2a..16596af885 100644 --- a/pvlib/tools.py +++ b/pvlib/tools.py @@ -2,8 +2,9 @@ Collection of functions used in pvlib_python """ +from collections import namedtuple import datetime as dt - +import warnings import numpy as np import pandas as pd import pytz @@ -251,3 +252,179 @@ def _build_kwargs(keys, input_dict): pass return kwargs + + +# FIXME: remove _array_newton when SciPy-1.2.0 is released +# pvlib.singlediode_methods.bishop88_i_from_v(..., method='newton') and other +# functions in singlediode_methods call scipy.optimize.newton with a vector +# unfortunately wrapping the functions with np.vectorize() was too slow +# a vectorized newton method was merged into SciPy but isn't released yet, so +# in the meantime, we just copied the relevant code: "_array_newton" for more +# info see: https://github.com/scipy/scipy/pull/8357 + +def _array_newton(func, x0, fprime, args, tol, maxiter, fprime2, + converged=False): + """ + A vectorized version of Newton, Halley, and secant methods for arrays. Do + not use this method directly. This method is called from :func:`newton` + when ``np.isscalar(x0)`` is true. For docstring, see :func:`newton`. + """ + try: + p = np.asarray(x0, dtype=float) + except TypeError: # can't convert complex to float + p = np.asarray(x0) + failures = np.ones_like(p, dtype=bool) # at start, nothing converged + nz_der = np.copy(failures) + if fprime is not None: + # Newton-Raphson method + for iteration in range(maxiter): + # first evaluate fval + fval = np.asarray(func(p, *args)) + # If all fval are 0, all roots have been found, then terminate + if not fval.any(): + failures = fval.astype(bool) + break + fder = np.asarray(fprime(p, *args)) + nz_der = (fder != 0) + # stop iterating if all derivatives are zero + if not nz_der.any(): + break + # Newton step + dp = fval[nz_der] / fder[nz_der] + if fprime2 is not None: + fder2 = np.asarray(fprime2(p, *args)) + dp = dp / (1.0 - 0.5 * dp * fder2[nz_der] / fder[nz_der]) + # only update nonzero derivatives + p[nz_der] -= dp + failures[nz_der] = np.abs(dp) >= tol # items not yet converged + # stop iterating if there aren't any failures, not incl zero der + if not failures[nz_der].any(): + break + else: + # Secant method + dx = np.finfo(float).eps**0.33 + p1 = p * (1 + dx) + np.where(p >= 0, dx, -dx) + q0 = np.asarray(func(p, *args)) + q1 = np.asarray(func(p1, *args)) + active = np.ones_like(p, dtype=bool) + for iteration in range(maxiter): + nz_der = (q1 != q0) + # stop iterating if all derivatives are zero + if not nz_der.any(): + p = (p1 + p) / 2.0 + break + # Secant Step + dp = (q1 * (p1 - p))[nz_der] / (q1 - q0)[nz_der] + # only update nonzero derivatives + p[nz_der] = p1[nz_der] - dp + active_zero_der = ~nz_der & active + p[active_zero_der] = (p1 + p)[active_zero_der] / 2.0 + active &= nz_der # don't assign zero derivatives again + failures[nz_der] = np.abs(dp) >= tol # not yet converged + # stop iterating if there aren't any failures, not incl zero der + if not failures[nz_der].any(): + break + p1, p = p, p1 + q0 = q1 + q1 = np.asarray(func(p1, *args)) + zero_der = ~nz_der & failures # don't include converged with zero-ders + if zero_der.any(): + # secant warnings + if fprime is None: + nonzero_dp = (p1 != p) + # non-zero dp, but infinite newton step + zero_der_nz_dp = (zero_der & nonzero_dp) + if zero_der_nz_dp.any(): + rms = np.sqrt( + sum((p1[zero_der_nz_dp] - p[zero_der_nz_dp]) ** 2) + ) + warnings.warn('RMS of {:g} reached'.format(rms), RuntimeWarning) + # newton or halley warnings + else: + all_or_some = 'all' if zero_der.all() else 'some' + msg = '{:s} derivatives were zero'.format(all_or_some) + warnings.warn(msg, RuntimeWarning) + elif failures.any(): + all_or_some = 'all' if failures.all() else 'some' + msg = '{0:s} failed to converge after {1:d} iterations'.format( + all_or_some, maxiter + ) + if failures.all(): + raise RuntimeError(msg) + warnings.warn(msg, RuntimeWarning) + if converged: + result = namedtuple('result', ('root', 'converged', 'zero_der')) + p = result(p, ~failures, zero_der) + return p + + +# Created April,2014 +# Author: Rob Andrews, Calama Consulting + +def _golden_sect_DataFrame(params, VL, VH, func): + """ + Vectorized golden section search for finding MPP from a dataframe + timeseries. + + Parameters + ---------- + params : dict + Dictionary containing scalars or arrays + of inputs to the function to be optimized. + Each row should represent an independent optimization. + + VL: float + Lower bound of the optimization + + VH: float + Upper bound of the optimization + + func: function + Function to be optimized must be in the form f(array-like, x) + + Returns + ------- + func(df,'V1') : DataFrame + function evaluated at the optimal point + + df['V1']: Dataframe + Dataframe of optimal points + + Notes + ----- + This function will find the MAXIMUM of a function + """ + + df = params + df['VH'] = VH + df['VL'] = VL + + err = df['VH'] - df['VL'] + errflag = True + iterations = 0 + + while errflag: + + phi = (np.sqrt(5)-1)/2*(df['VH']-df['VL']) + df['V1'] = df['VL'] + phi + df['V2'] = df['VH'] - phi + + df['f1'] = func(df, 'V1') + df['f2'] = func(df, 'V2') + df['SW_Flag'] = df['f1'] > df['f2'] + + df['VL'] = df['V2']*df['SW_Flag'] + df['VL']*(~df['SW_Flag']) + df['VH'] = df['V1']*~df['SW_Flag'] + df['VH']*(df['SW_Flag']) + + err = df['V1'] - df['V2'] + try: + errflag = (abs(err) > .01).any() + except ValueError: + errflag = (abs(err) > .01) + + iterations += 1 + + if iterations > 50: + raise Exception("EXCEPTION:iterations exceeded maximum (50)") + + return func(df, 'V1'), df['V1']