0
# Loss Models
1

2
Model various loss mechanisms that affect photovoltaic system performance including soiling, snow coverage, shading effects, and system losses. Comprehensive tools for quantifying and predicting energy losses under real-world conditions.
3

4
## Capabilities
5

6
### Soiling Models
7

8
Model the accumulation and cleaning of dust and debris on PV modules.
9

10
```python { .api }
11
def hsu(rainfall, cleaning_threshold, surface_tilt, pm2_5, pm10, 
12
        depo_veloc=None, rain_accum_period=None):
13
    """
14
    Calculate soiling losses using Hsu soiling model.
15
    
16
    Parameters:
17
    - rainfall: array-like, daily rainfall in mm
18
    - cleaning_threshold: numeric, rainfall threshold for cleaning in mm  
19
    - surface_tilt: numeric, surface tilt angle in degrees
20
    - pm2_5: array-like, PM2.5 concentration in μg/m³
21
    - pm10: array-like, PM10 concentration in μg/m³
22
    - depo_veloc: dict, deposition velocities by particle size
23
    - rain_accum_period: int, period for rainfall accumulation in days
24
    
25
    Returns:
26
    pandas.Series with soiling ratio (0-1, where 1 = clean)
27
    """
28

29
def kimber(rainfall, cleaning_threshold=6, soiling_loss_rate=0.0015,
30
           grace_period=14, max_loss_rate=0.3, manual_wash_dates=None):
31
    """
32
    Calculate soiling losses using Kimber soiling model.
33
    
34
    Parameters:
35
    - rainfall: pandas.Series, daily rainfall in mm
36
    - cleaning_threshold: numeric, rainfall threshold for cleaning in mm
37
    - soiling_loss_rate: numeric, daily soiling loss rate (fraction/day)
38
    - grace_period: int, days before soiling begins after cleaning
39
    - max_loss_rate: numeric, maximum soiling loss (0-1)
40
    - manual_wash_dates: array-like, dates of manual cleaning
41
    
42
    Returns:
43
    pandas.Series with soiling ratio (0-1, where 1 = clean)
44
    """
45
```
46

47
### Snow Coverage Models
48

49
Model snow accumulation and melting effects on PV modules.
50

51
```python { .api }
52
def fully_covered_nrel(snowfall, snow_depth=None, threshold_snowfall=1.0,
53
                       threshold_snow_depth=None):
54
    """
55
    Calculate snow coverage using NREL fully covered model.
56
    
57
    Parameters:
58
    - snowfall: array-like, daily snowfall in cm
59
    - snow_depth: array-like, snow depth on ground in cm (optional)
60
    - threshold_snowfall: numeric, snowfall threshold for full coverage in cm
61
    - threshold_snow_depth: numeric, snow depth threshold in cm
62
    
63
    Returns:
64
    pandas.Series with snow coverage fraction (0-1)
65
    """
66

67
def coverage_nrel(snowfall, poa_irradiance, temp_air, surface_tilt,
68
                  threshold_snowfall=1.0, threshold_temp=2.0):
69
    """
70
    Calculate snow coverage using NREL coverage model with melting.
71
    
72
    Parameters:
73
    - snowfall: array-like, daily snowfall in cm
74
    - poa_irradiance: array-like, plane-of-array irradiance in W/m²
75
    - temp_air: array-like, air temperature in degrees C
76
    - surface_tilt: numeric, surface tilt angle in degrees
77
    - threshold_snowfall: numeric, snowfall threshold in cm
78
    - threshold_temp: numeric, temperature threshold for melting in °C
79
    
80
    Returns:
81
    pandas.Series with snow coverage fraction (0-1)
82
    """
83

84
def dc_loss_nrel(snow_coverage, num_strings):
85
    """
86
    Calculate DC power loss from snow coverage using NREL model.
87
    
88
    Parameters:
89
    - snow_coverage: array-like, snow coverage fraction (0-1)
90
    - num_strings: int, number of strings in the array
91
    
92
    Returns:
93
    pandas.Series with DC loss factor (0-1, where 0 = no loss)
94
    """
95

96
def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
97
                  wind_speed, temp_air=None):
98
    """
99
    Calculate snow losses using Townsend model.
100
    
101
    Parameters:
102
    - snow_total: numeric, total seasonal snowfall in cm
103
    - snow_events: int, number of snow events per season
104
    - surface_tilt: numeric, surface tilt angle in degrees
105
    - relative_humidity: numeric, average relative humidity (0-100)
106
    - wind_speed: numeric, average wind speed in m/s
107
    - temp_air: numeric, average air temperature in °C (optional)
108
    
109
    Returns:
110
    numeric, seasonal snow loss fraction (0-1)
111
    """
112
```
113

114
### Shading Models
115

116
Model shading effects from nearby objects and inter-row shading.
117

118
```python { .api }
119
def ground_angle(surface_tilt, gcr, slant_height):
120
    """
121
    Calculate ground coverage angle for shading analysis.
122
    
123
    Parameters:
124
    - surface_tilt: numeric, surface tilt angle in degrees
125
    - gcr: numeric, ground coverage ratio (0-1)
126
    - slant_height: numeric, slant height of module in meters
127
    
128
    Returns:
129
    numeric, ground coverage angle in degrees
130
    """
131

132
def masking_angle(surface_tilt, gcr, slant_height):
133
    """
134
    Calculate masking angle for row-to-row shading.
135
    
136
    Parameters:
137
    - surface_tilt: numeric, surface tilt angle in degrees  
138
    - gcr: numeric, ground coverage ratio (0-1)
139
    - slant_height: numeric, slant height of module in meters
140
    
141
    Returns:
142
    numeric, masking angle in degrees
143
    """
144

145
def masking_angle_passias(surface_tilt, gcr):
146
    """
147
    Calculate masking angle using Passias method.
148
    
149
    Parameters:
150
    - surface_tilt: numeric, surface tilt angle in degrees
151
    - gcr: numeric, ground coverage ratio (0-1)
152
    
153
    Returns:
154
    numeric, masking angle in degrees
155
    """
156

157
def sky_diffuse_passias(masking_angle):
158
    """
159
    Calculate sky diffuse fraction using Passias method.
160
    
161
    Parameters:
162
    - masking_angle: numeric, masking angle in degrees
163
    
164
    Returns:
165
    numeric, sky diffuse fraction (0-1)
166
    """
167

168
def projected_solar_zenith_angle(solar_zenith, solar_azimuth, axis_azimuth):
169
    """
170
    Calculate projected solar zenith angle for tracking systems.
171
    
172
    Parameters:
173
    - solar_zenith: numeric, solar zenith angle in degrees
174
    - solar_azimuth: numeric, solar azimuth angle in degrees  
175
    - axis_azimuth: numeric, tracker axis azimuth in degrees
176
    
177
    Returns:
178
    numeric, projected solar zenith angle in degrees
179
    """
180

181
def shaded_fraction1d(solar_zenith, solar_azimuth, axis_azimuth,
182
                      gcr, cross_axis_tilt=0.0):
183
    """
184
    Calculate shaded fraction for 1D tracking systems.
185
    
186
    Parameters:
187
    - solar_zenith: numeric, solar zenith angle in degrees
188
    - solar_azimuth: numeric, solar azimuth angle in degrees
189
    - axis_azimuth: numeric, tracker axis azimuth in degrees
190
    - gcr: numeric, ground coverage ratio (0-1)
191
    - cross_axis_tilt: numeric, cross-axis tilt in degrees
192
    
193
    Returns:
194
    numeric, shaded fraction (0-1)
195
    """
196
```
197

198
### System Loss Factors
199

200
Model various electrical and optical system losses.
201

202
```python { .api }
203
def pvwatts_losses(soiling=2, shading=3, snow=0, mismatch=2, wiring=2,
204
                   connections=0.5, lid=1.5, nameplate_rating=1, age=0,
205
                   availability=3):
206
    """
207
    Calculate PVWatts system loss factors.
208
    
209
    Parameters:
210
    - soiling: numeric, soiling loss percentage
211
    - shading: numeric, shading loss percentage  
212
    - snow: numeric, snow loss percentage
213
    - mismatch: numeric, module mismatch loss percentage
214
    - wiring: numeric, DC wiring loss percentage
215
    - connections: numeric, DC connections loss percentage
216
    - lid: numeric, light-induced degradation percentage
217
    - nameplate_rating: numeric, nameplate rating loss percentage
218
    - age: numeric, age-related degradation percentage
219
    - availability: numeric, system availability loss percentage
220
    
221
    Returns:
222
    dict with individual loss factors and total system derate factor
223
    """
224

225
def dc_ohms_from_percent(vmp_ref, imp_ref, dc_ohmic_percent, 
226
                         modules_per_string=1, strings=1):
227
    """
228
    Convert DC loss percentage to resistance value.
229
    
230
    Parameters:
231
    - vmp_ref: numeric, module voltage at max power in volts
232
    - imp_ref: numeric, module current at max power in amps  
233
    - dc_ohmic_percent: numeric, DC ohmic loss percentage
234
    - modules_per_string: int, number of modules per string
235
    - strings: int, number of parallel strings
236
    
237
    Returns:
238
    numeric, equivalent series resistance in ohms
239
    """
240

241
def dc_ohmic_losses(resistance, current):
242
    """
243
    Calculate DC ohmic losses from resistance and current.
244
    
245
    Parameters:
246
    - resistance: numeric, series resistance in ohms
247
    - current: array-like, DC current in amps
248
    
249
    Returns:
250
    array-like, power loss in watts
251
    """
252

253
def combine_loss_factors(index, *losses, fill_method='ffill'):
254
    """
255
    Combine multiple loss factors into single time series.
256
    
257
    Parameters:
258
    - index: pandas.Index, time index for output
259
    - losses: multiple pandas.Series or scalars, loss factors to combine
260
    - fill_method: str, method to fill missing values
261
    
262
    Returns:
263
    pandas.Series with combined loss factors
264
    """
265
```
266

267
### Advanced Loss Models
268

269
Additional loss models for specific conditions and applications.
270

271
```python { .api }
272
def thermal_loss_factor(cell_temperature, temp_ref=25, temp_coeff=-0.004):
273
    """
274
    Calculate thermal loss factor from cell temperature.
275
    
276
    Parameters:
277
    - cell_temperature: array-like, cell temperature in °C
278
    - temp_ref: numeric, reference temperature in °C
279
    - temp_coeff: numeric, temperature coefficient per °C
280
    
281
    Returns:
282
    array-like, thermal loss factor (multiplicative, <1 indicates loss)
283
    """
284

285
def spectral_loss_factor(airmass, precipitable_water, module_type='csi'):
286
    """
287
    Calculate spectral mismatch loss factor.
288
    
289
    Parameters:
290
    - airmass: array-like, absolute airmass  
291
    - precipitable_water: array-like, precipitable water in cm
292
    - module_type: str, module technology type
293
    
294
    Returns:
295
    array-like, spectral loss factor (multiplicative)
296
    """
297

298
def iam_loss_factor(aoi, iam_model='physical', **kwargs):
299
    """
300
    Calculate incidence angle modifier loss factor.
301
    
302
    Parameters:
303
    - aoi: array-like, angle of incidence in degrees
304
    - iam_model: str, IAM model to use
305
    - kwargs: additional parameters for IAM model
306
    
307
    Returns:
308
    array-like, IAM loss factor (multiplicative, <1 indicates loss)
309
    """
310

311
def degradation_loss(years_operation, degradation_rate=0.005):
312
    """
313
    Calculate cumulative degradation losses over time.
314
    
315
    Parameters:
316
    - years_operation: numeric, years since installation
317
    - degradation_rate: numeric, annual degradation rate (fraction/year)
318
    
319
    Returns:
320
    numeric, cumulative degradation factor (<1 indicates loss)
321
    """
322
```
323

324
## Usage Examples
325

326
### Comprehensive Soiling Analysis
327

328
```python
329
import pvlib
330
from pvlib import soiling, iotools
331
import pandas as pd
332
import numpy as np
333
import matplotlib.pyplot as plt
334

335
# Load weather data with precipitation
336
lat, lon = 37.8756, -122.2441  # Berkeley, CA
337
try:
338
    # Try to get actual weather data
339
    weather_data, metadata = iotools.get_psm3(
340
        lat, lon, api_key='DEMO_KEY', email='user@example.com',
341
        names=2023, attributes=['ghi', 'temp_air', 'precip']
342
    )
343
    has_precip = 'precip' in weather_data.columns
344
except:
345
    # Generate synthetic weather data
346
    print("Using synthetic weather data")
347
    dates = pd.date_range('2023-01-01', '2023-12-31', freq='D')
348
    
349
    # Simulate seasonal precipitation pattern (more in winter)
350
    np.random.seed(42)
351
    seasonal_pattern = 1 + 0.8 * np.sin(2 * np.pi * (dates.dayofyear + 90) / 365)
352
    base_precip = np.random.exponential(2.0, len(dates))  # Exponential distribution
353
    daily_precip = base_precip * seasonal_pattern
354
    daily_precip[np.random.random(len(dates)) > 0.3] = 0  # 70% dry days
355
    
356
    weather_data = pd.DataFrame({
357
        'ghi': 200 + 300 * np.sin(2 * np.pi * dates.dayofyear / 365),  # Seasonal GHI
358
        'temp_air': 15 + 10 * np.sin(2 * np.pi * dates.dayofyear / 365),  # Seasonal temperature
359
        'precip': daily_precip
360
    }, index=dates)
361
    
362
    has_precip = True
363

364
if not has_precip:
365
    # Add synthetic precipitation if not available
366
    dates = weather_data.index
367
    np.random.seed(42)
368
    seasonal_pattern = 1 + 0.8 * np.sin(2 * np.pi * (dates.dayofyear + 90) / 365)
369
    base_precip = np.random.exponential(2.0, len(dates))
370
    daily_precip = base_precip * seasonal_pattern
371
    daily_precip[np.random.random(len(dates)) > 0.3] = 0
372
    weather_data['precip'] = daily_precip
373

374
# System parameters
375
surface_tilt = 25  # degrees
376
cleaning_threshold_low = 3   # mm (frequent cleaning)
377
cleaning_threshold_high = 10  # mm (infrequent cleaning)
378

379
# Soiling model parameters for different environments
380
environments = {
381
    'Clean': {'loss_rate': 0.0005, 'pm2_5': 8, 'pm10': 15},      # Rural/coastal
382
    'Moderate': {'loss_rate': 0.0015, 'pm2_5': 15, 'pm10': 30},  # Suburban
383
    'Dusty': {'loss_rate': 0.003, 'pm2_5': 30, 'pm10': 60}       # Desert/industrial
384
}
385

386
# Calculate soiling for different scenarios
387
soiling_results = {}
388

389
for env_name, params in environments.items():
390
    print(f"Calculating soiling for {env_name} environment...")
391
    
392
    # Kimber model with different cleaning thresholds
393
    for threshold in [cleaning_threshold_low, cleaning_threshold_high]:
394
        scenario_name = f"{env_name}_T{threshold}"
395
        
396
        soiling_ratio = soiling.kimber(
397
            rainfall=weather_data['precip'],
398
            cleaning_threshold=threshold,
399
            soiling_loss_rate=params['loss_rate'],
400
            grace_period=7,
401
            max_loss_rate=0.25
402
        )
403
        
404
        soiling_results[scenario_name] = soiling_ratio
405

406
    # Hsu model (if PM data available)
407
    if 'pm2_5' in params and 'pm10' in params:
408
        try:
409
            # Create constant PM concentrations (in real application, use time series)
410
            pm2_5_series = pd.Series(params['pm2_5'], index=weather_data.index)
411
            pm10_series = pd.Series(params['pm10'], index=weather_data.index)
412
            
413
            hsu_soiling = soiling.hsu(
414
                rainfall=weather_data['precip'],
415
                cleaning_threshold=cleaning_threshold_low,
416
                surface_tilt=surface_tilt,
417
                pm2_5=pm2_5_series,
418
                pm10=pm10_series
419
            )
420
            
421
            soiling_results[f"{env_name}_Hsu"] = hsu_soiling
422
        except Exception as e:
423
            print(f"Hsu model not available: {e}")
424

425
# Calculate energy impact
426
poa_clean = weather_data['ghi'] * 1.2  # Simplified POA calculation
427
annual_energy_impacts = {}
428

429
for scenario, soiling_ratio in soiling_results.items():
430
    # Apply soiling losses
431
    poa_soiled = poa_clean * soiling_ratio
432
    
433
    # Calculate energy totals
434
    clean_energy = poa_clean.sum() / 1000  # kWh/m²/year
435
    soiled_energy = poa_soiled.sum() / 1000
436
    
437
    energy_loss = (clean_energy - soiled_energy) / clean_energy * 100
438
    
439
    annual_energy_impacts[scenario] = {
440
        'clean_energy': clean_energy,
441
        'soiled_energy': soiled_energy, 
442
        'energy_loss_pct': energy_loss,
443
        'avg_soiling_ratio': soiling_ratio.mean()
444
    }
445

446
# Plot soiling analysis
447
fig, axes = plt.subplots(3, 2, figsize=(15, 12))
448

449
# Time series comparison for one environment
450
env_scenarios = [k for k in soiling_results.keys() if k.startswith('Moderate')]
451
for scenario in env_scenarios:
452
    axes[0, 0].plot(soiling_results[scenario].index, 
453
                   soiling_results[scenario] * 100,
454
                   label=scenario.replace('Moderate_', ''), linewidth=1.5)
455

456
axes[0, 0].set_ylabel('Soiling Ratio (%)')
457
axes[0, 0].set_title('Soiling Time Series - Moderate Environment')
458
axes[0, 0].legend()
459
axes[0, 0].grid(True)
460

461
# Monthly averages comparison
462
monthly_soiling = {}
463
for scenario in soiling_results:
464
    monthly_avg = soiling_results[scenario].resample('M').mean() * 100
465
    monthly_soiling[scenario] = monthly_avg
466

467
# Plot monthly averages for clean vs dusty environments
468
clean_scenarios = [k for k in monthly_soiling.keys() if 'Clean' in k]
469
dusty_scenarios = [k for k in monthly_soiling.keys() if 'Dusty' in k]
470

471
for scenario in clean_scenarios:
472
    axes[0, 1].plot(monthly_soiling[scenario].index.month, 
473
                   monthly_soiling[scenario], 
474
                   label=scenario, linewidth=2, linestyle='-')
475

476
for scenario in dusty_scenarios:
477
    axes[0, 1].plot(monthly_soiling[scenario].index.month, 
478
                   monthly_soiling[scenario], 
479
                   label=scenario, linewidth=2, linestyle='--')
480

481
axes[0, 1].set_xlabel('Month')
482
axes[0, 1].set_ylabel('Avg Soiling Ratio (%)')
483
axes[0, 1].set_title('Monthly Average Soiling Ratios')
484
axes[0, 1].legend()
485
axes[0, 1].grid(True)
486

487
# Energy impact summary
488
scenarios = list(annual_energy_impacts.keys())
489
energy_losses = [annual_energy_impacts[s]['energy_loss_pct'] for s in scenarios]
490
avg_soiling = [annual_energy_impacts[s]['avg_soiling_ratio'] * 100 for s in scenarios]
491

492
colors = ['blue' if 'Clean' in s else 'green' if 'Moderate' in s else 'red' 
493
         for s in scenarios]
494

495
bars = axes[1, 0].bar(range(len(scenarios)), energy_losses, color=colors, alpha=0.7)
496
axes[1, 0].set_xticks(range(len(scenarios)))
497
axes[1, 0].set_xticklabels(scenarios, rotation=45, ha='right')
498
axes[1, 0].set_ylabel('Annual Energy Loss (%)')
499
axes[1, 0].set_title('Annual Energy Impact by Scenario')
500
axes[1, 0].grid(True, axis='y')
501

502
# Add value labels on bars
503
for bar, loss in zip(bars, energy_losses):
504
    height = bar.get_height()
505
    axes[1, 0].text(bar.get_x() + bar.get_width()/2., height + 0.1,
506
                   f'{loss:.1f}%', ha='center', va='bottom')
507

508
# Correlation: precipitation vs soiling recovery
509
precip_monthly = weather_data['precip'].resample('M').sum()
510
# Use moderate environment with low threshold as reference
511
ref_soiling = monthly_soiling['Moderate_T3']
512

513
axes[1, 1].scatter(precip_monthly, ref_soiling, alpha=0.7, s=50)
514
axes[1, 1].set_xlabel('Monthly Precipitation (mm)')
515
axes[1, 1].set_ylabel('Monthly Avg Soiling Ratio (%)')
516
axes[1, 1].set_title('Precipitation vs Soiling Recovery')
517
axes[1, 1].grid(True)
518

519
# Add trend line
520
z = np.polyfit(precip_monthly, ref_soiling, 1)
521
p = np.poly1d(z)
522
axes[1, 1].plot(precip_monthly, p(precip_monthly), "r--", alpha=0.8, linewidth=2)
523

524
# Cleaning event analysis
525
cleaning_events = (weather_data['precip'] >= cleaning_threshold_low).sum()
526
days_between_cleaning = 365 / cleaning_events if cleaning_events > 0 else 365
527

528
# Economic analysis (simplified)
529
system_size = 100  # kW
530
electricity_rate = 0.15  # $/kWh
531
cleaning_cost_per_event = 200  # $
532

533
economic_analysis = []
534
for scenario, results in annual_energy_impacts.items():
535
    threshold = int(scenario.split('_T')[-1]) if '_T' in scenario else cleaning_threshold_low
536
    
537
    # Annual cleaning events
538
    annual_events = (weather_data['precip'] >= threshold).sum()
539
    
540
    # Annual costs and savings
541
    energy_loss_kwh = results['energy_loss_pct'] / 100 * system_size * 1500  # Assume 1500 kWh/kW/year
542
    revenue_loss = energy_loss_kwh * electricity_rate
543
    cleaning_costs = annual_events * cleaning_cost_per_event
544
    
545
    economic_analysis.append({
546
        'scenario': scenario,
547
        'threshold': threshold,
548
        'events': annual_events,
549
        'energy_loss_pct': results['energy_loss_pct'],
550
        'revenue_loss': revenue_loss,
551
        'cleaning_costs': cleaning_costs,
552
        'net_cost': revenue_loss + cleaning_costs
553
    })
554

555
# Plot economic analysis
556
econ_df = pd.DataFrame(economic_analysis)
557
econ_summary = econ_df.groupby('threshold').agg({
558
    'energy_loss_pct': 'mean',
559
    'revenue_loss': 'mean', 
560
    'cleaning_costs': 'mean',
561
    'net_cost': 'mean'
562
}).reset_index()
563

564
x_pos = np.arange(len(econ_summary))
565
width = 0.35
566

567
bars1 = axes[2, 0].bar(x_pos - width/2, econ_summary['revenue_loss'], 
568
                      width, label='Revenue Loss', alpha=0.7, color='red')
569
bars2 = axes[2, 0].bar(x_pos + width/2, econ_summary['cleaning_costs'], 
570
                      width, label='Cleaning Costs', alpha=0.7, color='blue')
571

572
axes[2, 0].set_xlabel('Cleaning Threshold (mm)')
573
axes[2, 0].set_ylabel('Annual Cost ($)')
574
axes[2, 0].set_title('Economic Impact of Cleaning Strategy')
575
axes[2, 0].set_xticks(x_pos)
576
axes[2, 0].set_xticklabels(econ_summary['threshold'])
577
axes[2, 0].legend()
578
axes[2, 0].grid(True, axis='y')
579

580
# Optimal cleaning threshold
581
axes[2, 1].plot(econ_summary['threshold'], econ_summary['net_cost'], 'go-', 
582
               linewidth=3, markersize=8)
583
axes[2, 1].set_xlabel('Cleaning Threshold (mm)')
584
axes[2, 1].set_ylabel('Total Annual Cost ($)')
585
axes[2, 1].set_title('Total Cost vs Cleaning Threshold')
586
axes[2, 1].grid(True)
587

588
# Find and mark optimal threshold
589
optimal_idx = econ_summary['net_cost'].idxmin()
590
optimal_threshold = econ_summary.loc[optimal_idx, 'threshold']
591
optimal_cost = econ_summary.loc[optimal_idx, 'net_cost']
592

593
axes[2, 1].scatter([optimal_threshold], [optimal_cost], 
594
                  color='red', s=100, zorder=5)
595
axes[2, 1].annotate(f'Optimal: {optimal_threshold}mm\n${optimal_cost:.0f}/year',
596
                   xy=(optimal_threshold, optimal_cost),
597
                   xytext=(optimal_threshold + 1, optimal_cost + 200),
598
                   arrowprops=dict(arrowstyle='->', color='red'),
599
                   fontsize=10, fontweight='bold')
600

601
plt.tight_layout()
602
plt.show()
603

604
# Print comprehensive results
605
print(f"\nComprehensive Soiling Analysis Results")
606
print("="*50)
607
print(f"Location: {lat:.2f}°N, {lon:.2f}°W")
608
print(f"Annual precipitation: {weather_data['precip'].sum():.1f} mm")
609
print(f"Cleaning events (3mm threshold): {(weather_data['precip'] >= 3).sum()}")
610
print(f"Cleaning events (10mm threshold): {(weather_data['precip'] >= 10).sum()}")
611

612
print(f"\nEnergy Loss Summary:")
613
for scenario, results in annual_energy_impacts.items():
614
    print(f"{scenario:15s}: {results['energy_loss_pct']:5.2f}% "
615
          f"(avg soiling ratio: {results['avg_soiling_ratio']*100:.1f}%)")
616

617
print(f"\nOptimal Cleaning Strategy:")
618
print(f"Threshold: {optimal_threshold} mm")
619
print(f"Total annual cost: ${optimal_cost:.0f}")
620
print(f"Revenue loss: ${econ_summary.loc[optimal_idx, 'revenue_loss']:.0f}")
621
print(f"Cleaning costs: ${econ_summary.loc[optimal_idx, 'cleaning_costs']:.0f}")
622
```
623

624
### Snow Loss Analysis
625

626
```python
627
import pvlib
628
from pvlib import snow, iotools, solarposition, irradiance
629
import pandas as pd
630
import numpy as np
631
import matplotlib.pyplot as plt
632

633
# Location with significant snowfall
634
lat, lon = 44.0582, -121.3153  # Bend, Oregon
635
elevation = 1116  # meters
636

637
# Generate winter weather data (simplified)
638
winter_dates = pd.date_range('2022-11-01', '2023-03-31', freq='D')
639
np.random.seed(42)
640

641
# Simulate winter conditions
642
temp_base = -2 + 8 * np.sin(2 * np.pi * winter_dates.dayofyear / 365)
643
temp_noise = np.random.normal(0, 3, len(winter_dates))
644
temp_air = temp_base + temp_noise
645

646
# Snowfall occurs when temp < 2°C
647
snowfall_prob = 1 / (1 + np.exp(2 * (temp_air - 0)))  # Sigmoid probability
648
snowfall_events = np.random.random(len(winter_dates)) < snowfall_prob * 0.3
649
snowfall_amounts = np.where(snowfall_events, 
650
                           np.random.exponential(3, len(winter_dates)), 0)
651

652
# Snow depth simulation (simplified accumulation/melt)
653
snow_depth = np.zeros(len(winter_dates))
654
for i in range(1, len(winter_dates)):
655
    # Accumulation
656
    snow_depth[i] = snow_depth[i-1] + snowfall_amounts[i]
657
    
658
    # Melting when temp > 0°C
659
    if temp_air[i] > 0:
660
        melt_rate = 0.3 * temp_air[i]  # cm/day per °C above 0
661
        snow_depth[i] = max(0, snow_depth[i] - melt_rate)
662

663
# Simulate solar irradiance
664
clear_sky_winter = pvlib.clearsky.ineichen(
665
    winter_dates, lat, lon, altitude=elevation
666
)
667

668
# System parameters
669
surface_tilt = 35  # degrees (steeper for snow shedding)
670
surface_azimuth = 180  # south-facing
671

672
# Calculate solar position and POA irradiance
673
solar_pos = solarposition.get_solarposition(winter_dates, lat, lon)
674
poa_irradiance = irradiance.get_total_irradiance(
675
    surface_tilt, surface_azimuth,
676
    solar_pos['zenith'], solar_pos['azimuth'],
677
    clear_sky_winter['dni'], clear_sky_winter['ghi'], clear_sky_winter['dhi']
678
)['poa_global']
679

680
# Create weather DataFrame
681
weather_winter = pd.DataFrame({
682
    'temp_air': temp_air,
683
    'snowfall': snowfall_amounts, 
684
    'snow_depth': snow_depth,
685
    'poa_irradiance': poa_irradiance
686
}, index=winter_dates)
687

688
# Snow models comparison
689
print("Calculating snow coverage models...")
690

691
# NREL fully covered model
692
snow_coverage_basic = snow.fully_covered_nrel(
693
    snowfall=weather_winter['snowfall'],
694
    threshold_snowfall=1.0  # 1 cm threshold
695
)
696

697
# NREL coverage model with melting
698
snow_coverage_advanced = snow.coverage_nrel(
699
    snowfall=weather_winter['snowfall'],
700
    poa_irradiance=weather_winter['poa_irradiance'],
701
    temp_air=weather_winter['temp_air'],
702
    surface_tilt=surface_tilt,
703
    threshold_snowfall=1.0,
704
    threshold_temp=2.0
705
)
706

707
# Townsend seasonal model
708
seasonal_snow_total = weather_winter['snowfall'].sum()
709
seasonal_snow_events = (weather_winter['snowfall'] > 0).sum()
710

711
townsend_loss = snow.loss_townsend(
712
    snow_total=seasonal_snow_total,
713
    snow_events=seasonal_snow_events,
714
    surface_tilt=surface_tilt,
715
    relative_humidity=70,  # Assumed average
716
    wind_speed=3.0  # Assumed average m/s
717
)
718

719
print(f"Townsend seasonal snow loss: {townsend_loss*100:.1f}%")
720

721
# Calculate energy impacts
722
num_strings = 10  # Example system configuration
723

724
# DC losses from snow coverage
725
dc_loss_basic = snow.dc_loss_nrel(snow_coverage_basic, num_strings)
726
dc_loss_advanced = snow.dc_loss_nrel(snow_coverage_advanced, num_strings)
727

728
# Energy calculations
729
baseline_energy = weather_winter['poa_irradiance'].sum() / 1000  # kWh/m²
730
energy_basic = (weather_winter['poa_irradiance'] * (1 - dc_loss_basic)).sum() / 1000
731
energy_advanced = (weather_winter['poa_irradiance'] * (1 - dc_loss_advanced)).sum() / 1000
732

733
energy_loss_basic = (baseline_energy - energy_basic) / baseline_energy * 100
734
energy_loss_advanced = (baseline_energy - energy_advanced) / baseline_energy * 100
735

736
print(f"\nWinter Energy Analysis:")
737
print(f"Baseline energy (no snow): {baseline_energy:.1f} kWh/m²")
738
print(f"Basic model energy: {energy_basic:.1f} kWh/m² (loss: {energy_loss_basic:.1f}%)")
739
print(f"Advanced model energy: {energy_advanced:.1f} kWh/m² (loss: {energy_loss_advanced:.1f}%)")
740

741
# Plot snow analysis
742
fig, axes = plt.subplots(4, 1, figsize=(15, 12))
743

744
# Weather conditions
745
axes[0].plot(weather_winter.index, weather_winter['temp_air'], 
746
            'b-', label='Air Temperature', linewidth=1.5)
747
axes[0].axhline(y=0, color='k', linestyle='--', alpha=0.5)
748
axes[0].set_ylabel('Temperature (°C)')
749
axes[0].set_title('Winter Weather Conditions')
750
axes[0].legend()
751
axes[0].grid(True)
752

753
# Add snowfall as bars on secondary axis
754
ax0_twin = axes[0].twinx()
755
snow_events_idx = weather_winter['snowfall'] > 0
756
ax0_twin.bar(weather_winter.index[snow_events_idx], 
757
            weather_winter['snowfall'][snow_events_idx],
758
            width=1, alpha=0.3, color='lightblue', label='Snowfall')
759
ax0_twin.set_ylabel('Snowfall (cm)', color='lightblue')
760
ax0_twin.legend(loc='upper right')
761

762
# Snow depth and coverage
763
axes[1].plot(weather_winter.index, weather_winter['snow_depth'], 
764
            'g-', label='Snow Depth', linewidth=2)
765
axes[1].set_ylabel('Snow Depth (cm)')
766
axes[1].set_title('Snow Accumulation and Melt')
767
axes[1].legend()
768
axes[1].grid(True)
769

770
# Snow coverage comparison
771
axes[2].plot(weather_winter.index, snow_coverage_basic * 100, 
772
            'r-', label='Basic Model', linewidth=2)
773
axes[2].plot(weather_winter.index, snow_coverage_advanced * 100, 
774
            'b-', label='Advanced Model', linewidth=2)
775
axes[2].set_ylabel('Snow Coverage (%)')
776
axes[2].set_title('Snow Coverage Models Comparison')
777
axes[2].legend()
778
axes[2].grid(True)
779

780
# Energy impact
781
daily_energy_baseline = weather_winter['poa_irradiance'] / 1000
782
daily_energy_basic = weather_winter['poa_irradiance'] * (1 - dc_loss_basic) / 1000
783
daily_energy_advanced = weather_winter['poa_irradiance'] * (1 - dc_loss_advanced) / 1000
784

785
axes[3].plot(weather_winter.index, daily_energy_baseline, 
786
            'k--', label='No Snow', linewidth=2, alpha=0.7)
787
axes[3].plot(weather_winter.index, daily_energy_basic, 
788
            'r-', label='Basic Model', linewidth=1.5)
789
axes[3].plot(weather_winter.index, daily_energy_advanced, 
790
            'b-', label='Advanced Model', linewidth=1.5)
791
axes[3].set_ylabel('Daily Energy (kWh/m²)')
792
axes[3].set_xlabel('Date')
793
axes[3].set_title('Daily Energy Production with Snow Losses')
794
axes[3].legend()
795
axes[3].grid(True)
796

797
plt.tight_layout()
798
plt.show()
799

800
# Monthly analysis
801
monthly_stats = weather_winter.groupby(weather_winter.index.month).agg({
802
    'temp_air': ['mean', 'min', 'max'],
803
    'snowfall': 'sum',
804
    'snow_depth': 'mean',
805
    'poa_irradiance': 'sum'
806
}).round(2)
807

808
monthly_coverage = pd.DataFrame({
809
    'basic_coverage': snow_coverage_basic.groupby(snow_coverage_basic.index.month).mean(),
810
    'advanced_coverage': snow_coverage_advanced.groupby(snow_coverage_advanced.index.month).mean(),
811
    'basic_loss': dc_loss_basic.groupby(dc_loss_basic.index.month).mean(),
812
    'advanced_loss': dc_loss_advanced.groupby(dc_loss_advanced.index.month).mean()
813
}) * 100
814

815
print(f"\nMonthly Snow Impact Summary:")
816
print("="*40)
817
for month in monthly_coverage.index:
818
    month_name = pd.Timestamp(2023, month, 1).strftime('%B')
819
    print(f"{month_name}:")
820
    print(f"  Avg Coverage: Basic {monthly_coverage.loc[month, 'basic_coverage']:.1f}%, "
821
          f"Advanced {monthly_coverage.loc[month, 'advanced_coverage']:.1f}%")
822
    print(f"  Avg DC Loss:  Basic {monthly_coverage.loc[month, 'basic_loss']:.1f}%, "
823
          f"Advanced {monthly_coverage.loc[month, 'advanced_loss']:.1f}%")
824

825
# Tilt angle sensitivity analysis
826
tilt_angles = np.arange(15, 61, 5)
827
tilt_analysis = []
828

829
for tilt in tilt_angles:
830
    # Recalculate POA for this tilt
831
    poa_tilt = irradiance.get_total_irradiance(
832
        tilt, surface_azimuth,
833
        solar_pos['zenith'], solar_pos['azimuth'],
834
        clear_sky_winter['dni'], clear_sky_winter['ghi'], clear_sky_winter['dhi']
835
    )['poa_global']
836
    
837
    # Snow coverage with this tilt
838
    snow_coverage_tilt = snow.coverage_nrel(
839
        snowfall=weather_winter['snowfall'],
840
        poa_irradiance=poa_tilt,
841
        temp_air=weather_winter['temp_air'],
842
        surface_tilt=tilt,
843
        threshold_snowfall=1.0,
844
        threshold_temp=2.0
845
    )
846
    
847
    # Energy impact
848
    dc_loss_tilt = snow.dc_loss_nrel(snow_coverage_tilt, num_strings)
849
    energy_tilt = (poa_tilt * (1 - dc_loss_tilt)).sum() / 1000
850
    baseline_tilt = poa_tilt.sum() / 1000
851
    loss_pct = (baseline_tilt - energy_tilt) / baseline_tilt * 100
852
    
853
    tilt_analysis.append({
854
        'tilt': tilt,
855
        'baseline_energy': baseline_tilt,
856
        'actual_energy': energy_tilt,
857
        'snow_loss_pct': loss_pct,
858
        'avg_coverage': snow_coverage_tilt.mean() * 100
859
    })
860

861
tilt_df = pd.DataFrame(tilt_analysis)
862

863
print(f"\nTilt Angle Optimization for Snow:")
864
print("="*35)
865
optimal_tilt_idx = tilt_df['snow_loss_pct'].idxmin()
866
optimal_tilt = tilt_df.loc[optimal_tilt_idx, 'tilt']
867
min_snow_loss = tilt_df.loc[optimal_tilt_idx, 'snow_loss_pct']
868

869
print(f"Optimal tilt angle: {optimal_tilt}°")
870
print(f"Minimum snow loss: {min_snow_loss:.1f}%")
871

872
# Compare with latitude-based rule of thumb
873
latitude_tilt = lat
874
lat_idx = np.argmin(np.abs(tilt_df['tilt'] - latitude_tilt))
875
lat_snow_loss = tilt_df.loc[lat_idx, 'snow_loss_pct']
876

877
print(f"Latitude tilt ({latitude_tilt:.0f}°): {lat_snow_loss:.1f}% snow loss")
878
print(f"Improvement with optimal tilt: {lat_snow_loss - min_snow_loss:.1f} percentage points")
879
```
880

881
### Comprehensive Loss Factor Analysis
882

883
```python
884
import pvlib
885
from pvlib import pvsystem, temperature, atmosphere, solarposition
886
import pandas as pd
887
import numpy as np
888
import matplotlib.pyplot as plt
889

890
# System and location parameters
891
lat, lon = 40.0150, -105.2705  # Boulder, CO
892
system_params = {
893
    'module_power': 400,      # W
894
    'modules_per_string': 12,
895
    'strings': 8,
896
    'system_power': 38400,    # W total
897
    'surface_tilt': 30,       # degrees
898
    'surface_azimuth': 180,   # south
899
    'inverter_efficiency': 0.96
900
}
901

902
# Create annual time series
903
times = pd.date_range('2023-01-01', '2023-12-31 23:00', freq='H')
904

905
# Calculate solar position and clear sky
906
solar_pos = solarposition.get_solarposition(times, lat, lon)
907
clear_sky = pvlib.clearsky.ineichen(times, lat, lon, altitude=1655)
908

909
# Calculate POA irradiance
910
poa = pvlib.irradiance.get_total_irradiance(
911
    system_params['surface_tilt'], 
912
    system_params['surface_azimuth'],
913
    solar_pos['zenith'], 
914
    solar_pos['azimuth'],
915
    clear_sky['dni'], 
916
    clear_sky['ghi'], 
917
    clear_sky['dhi']
918
)
919

920
# Simulate weather variations
921
np.random.seed(42)
922
weather_variations = pd.DataFrame({
923
    'ghi_factor': 0.85 + 0.3 * np.random.random(len(times)),  # Cloud variations
924
    'temp_air': 15 + 20 * np.sin(2 * np.pi * times.dayofyear / 365) + np.random.normal(0, 3, len(times)),
925
    'wind_speed': 2 + 3 * np.random.exponential(1, len(times)),
926
    'relative_humidity': 30 + 40 * np.random.random(len(times)),
927
    'precipitation': np.random.exponential(1, len(times)) * (np.random.random(len(times)) < 0.1)
928
}, index=times)
929

930
# Apply weather variations to irradiance
931
poa_actual = poa['poa_global'] * weather_variations['ghi_factor']
932

933
# Calculate atmospheric parameters for spectral losses
934
airmass_rel = atmosphere.get_relative_airmass(solar_pos['zenith'])
935
airmass_abs = atmosphere.get_absolute_airmass(airmass_rel, pressure=85000)
936
precipitable_water = 1.5 + 0.5 * np.sin(2 * np.pi * times.dayofyear / 365)  # Seasonal variation
937

938
# Define all loss factors
939
print("Calculating comprehensive loss factors...")
940

941
# 1. PVWatts standard losses
942
pvwatts_loss_factors = pvsystem.pvwatts_losses(
943
    soiling=2.0,          # 2% soiling
944
    shading=1.0,          # 1% far shading  
945
    snow=0.5,             # 0.5% snow (minimal for this location)
946
    mismatch=2.0,         # 2% module mismatch
947
    wiring=2.0,           # 2% DC wiring
948
    connections=0.5,      # 0.5% connections
949
    lid=1.5,              # 1.5% light-induced degradation  
950
    nameplate_rating=1.0, # 1% nameplate rating
951
    age=0.0,              # New system
952
    availability=3.0      # 3% system availability
953
)
954

955
total_pvwatts_derate = pvwatts_loss_factors['total']
956
print(f"PVWatts total derate factor: {total_pvwatts_derate:.3f}")
957

958
# 2. Temperature losses
959
cell_temp = temperature.sapm_cell(
960
    poa_actual,
961
    weather_variations['temp_air'],
962
    weather_variations['wind_speed'],
963
    a=-3.47, b=-0.0594, deltaT=3  # Typical glass/cell/polymer values
964
)
965

966
temp_coeff = -0.004  # %/°C, typical for c-Si
967
temp_loss_factor = 1 + temp_coeff * (cell_temp - 25)
968

969
# 3. Spectral losses (simplified)
970
from pvlib import spectrum
971
sapm_spectral_params = {
972
    'a0': 0.928, 'a1': 0.068, 'a2': -0.0077, 
973
    'a3': 0.0001, 'a4': -0.000002
974
}
975
spectral_factor = spectrum.spectral_factor_sapm(airmass_abs, sapm_spectral_params)
976

977
# 4. Incidence angle losses
978
aoi = pvlib.irradiance.aoi(
979
    system_params['surface_tilt'],
980
    system_params['surface_azimuth'], 
981
    solar_pos['zenith'],
982
    solar_pos['azimuth']
983
)
984
iam_factor = pvlib.iam.ashrae(aoi, b=0.05)
985

986
# 5. DC ohmic losses
987
vmp_ref = 37.8  # V at STC
988
imp_ref = 10.58  # A at STC
989
dc_ohmic_percent = 1.5  # 1.5% loss
990

991
dc_resistance = pvsystem.dc_ohms_from_percent(
992
    vmp_ref, imp_ref, dc_ohmic_percent,
993
    modules_per_string=system_params['modules_per_string'],
994
    strings=system_params['strings']
995
)
996

997
# Simplified current calculation (would normally use full IV model)
998
estimated_current = poa_actual / 1000 * imp_ref * system_params['strings']
999
dc_ohmic_loss_factor = 1 - pvsystem.dc_ohmic_losses(dc_resistance, estimated_current) / (estimated_current * vmp_ref * system_params['modules_per_string'])
1000

1001
# 6. Soiling losses (time-varying)
1002
soiling_factor = pvsystem.soiling.kimber(
1003
    rainfall=weather_variations['precipitation'],
1004
    cleaning_threshold=5,  # mm
1005
    soiling_loss_rate=0.002,  # per day
1006
    max_loss_rate=0.15
1007
)
1008

1009
# 7. Degradation (age-related)
1010
system_age = 5  # years
1011
degradation_rate = 0.005  # 0.5% per year
1012
degradation_factor = (1 - degradation_rate) ** system_age
1013

1014
# Combine all loss factors
1015
combined_losses = pd.DataFrame({
1016
    'baseline_poa': poa_actual,
1017
    'temp_factor': temp_loss_factor,
1018
    'spectral_factor': spectral_factor,
1019
    'iam_factor': iam_factor, 
1020
    'dc_ohmic_factor': dc_ohmic_loss_factor,
1021
    'soiling_factor': soiling_factor,
1022
    'degradation_factor': degradation_factor,
1023
    'pvwatts_derate': total_pvwatts_derate
1024
}, index=times)
1025

1026
# Calculate cumulative effects
1027
combined_losses['effective_irradiance'] = (
1028
    combined_losses['baseline_poa'] *
1029
    combined_losses['temp_factor'] *
1030
    combined_losses['spectral_factor'] *
1031
    combined_losses['iam_factor'] *
1032
    combined_losses['dc_ohmic_factor'] *
1033
    combined_losses['soiling_factor'] *
1034
    combined_losses['degradation_factor'] *
1035
    combined_losses['pvwatts_derate']
1036
)
1037

1038
# Calculate individual loss impacts
1039
loss_impacts = {}
1040
for factor in ['temp_factor', 'spectral_factor', 'iam_factor', 
1041
               'dc_ohmic_factor', 'soiling_factor']:
1042
    annual_avg = combined_losses[factor].mean()
1043
    loss_impacts[factor] = (1 - annual_avg) * 100
1044

1045
# Add fixed losses
1046
loss_impacts['pvwatts_losses'] = (1 - total_pvwatts_derate) * 100
1047
loss_impacts['degradation'] = (1 - degradation_factor) * 100
1048

1049
# Energy analysis
1050
baseline_energy = combined_losses['baseline_poa'].sum() / 1000  # kWh/m²/year
1051
effective_energy = combined_losses['effective_irradiance'].sum() / 1000
1052
total_loss_pct = (baseline_energy - effective_energy) / baseline_energy * 100
1053

1054
print(f"\nAnnual Energy Analysis:")
1055
print("="*30)
1056
print(f"Baseline POA energy: {baseline_energy:.1f} kWh/m²/year")
1057
print(f"Effective energy: {effective_energy:.1f} kWh/m²/year") 
1058
print(f"Total energy loss: {total_loss_pct:.1f}%")
1059

1060
print(f"\nIndividual Loss Factor Impacts:")
1061
print("="*35)
1062
for factor, impact in sorted(loss_impacts.items(), key=lambda x: x[1], reverse=True):
1063
    factor_name = factor.replace('_factor', '').replace('_', ' ').title()
1064
    print(f"{factor_name:15s}: {impact:5.2f}%")
1065

1066
# Plot comprehensive loss analysis
1067
fig, axes = plt.subplots(4, 2, figsize=(16, 14))
1068

1069
# Time series of major loss factors
1070
sample_days = slice('2023-06-15', '2023-06-22')  # Summer week
1071
daily_data = combined_losses[sample_days]
1072

1073
axes[0, 0].plot(daily_data.index, daily_data['temp_factor'], 
1074
               label='Temperature', linewidth=2)
1075
axes[0, 0].plot(daily_data.index, daily_data['spectral_factor'], 
1076
               label='Spectral', linewidth=2)
1077
axes[0, 0].plot(daily_data.index, daily_data['iam_factor'], 
1078
               label='IAM', linewidth=2)
1079
axes[0, 0].set_ylabel('Loss Factor')
1080
axes[0, 0].set_title('Time-Varying Loss Factors (Summer Week)')
1081
axes[0, 0].legend()
1082
axes[0, 0].grid(True)
1083

1084
# Monthly averages
1085
monthly_factors = combined_losses.resample('M').mean()
1086
months = monthly_factors.index.month
1087

1088
axes[0, 1].plot(months, monthly_factors['temp_factor'], 'ro-', 
1089
               label='Temperature', linewidth=2)
1090
axes[0, 1].plot(months, monthly_factors['spectral_factor'], 'bs-', 
1091
               label='Spectral', linewidth=2)
1092
axes[0, 1].plot(months, monthly_factors['iam_factor'], 'g^-', 
1093
               label='IAM', linewidth=2)
1094
axes[0, 1].set_xlabel('Month')
1095
axes[0, 1].set_ylabel('Monthly Avg Factor')
1096
axes[0, 1].set_title('Seasonal Loss Factor Variations')
1097
axes[0, 1].legend()
1098
axes[0, 1].grid(True)
1099

1100
# Soiling time series
1101
axes[1, 0].plot(combined_losses.index, combined_losses['soiling_factor'], 
1102
               'brown', linewidth=1)
1103
axes[1, 0].set_ylabel('Soiling Factor')
1104
axes[1, 0].set_title('Soiling Losses Throughout Year')
1105
axes[1, 0].grid(True)
1106

1107
# Add precipitation events
1108
precip_events = weather_variations['precipitation'] > 3
1109
axes[1, 0].scatter(weather_variations.index[precip_events], 
1110
                  [0.95] * precip_events.sum(),
1111
                  c='blue', alpha=0.6, s=10, label='Rain Events')
1112
axes[1, 0].legend()
1113

1114
# Loss impact waterfall chart
1115
factors = list(loss_impacts.keys())
1116
impacts = list(loss_impacts.values())
1117
cumulative = np.cumsum([0] + impacts)
1118

1119
bars = axes[1, 1].bar(range(len(factors)), impacts, 
1120
                     color='lightcoral', alpha=0.7)
1121
axes[1, 1].set_xticks(range(len(factors)))
1122
axes[1, 1].set_xticklabels([f.replace('_', '\n') for f in factors], 
1123
                          rotation=45, ha='right')
1124
axes[1, 1].set_ylabel('Loss Impact (%)')
1125
axes[1, 1].set_title('Individual Loss Impacts')
1126
axes[1, 1].grid(True, axis='y')
1127

1128
# Add value labels
1129
for bar, impact in zip(bars, impacts):
1130
    height = bar.get_height()
1131
    axes[1, 1].text(bar.get_x() + bar.get_width()/2., height + 0.1,
1132
                   f'{impact:.1f}%', ha='center', va='bottom', fontsize=9)
1133

1134
# Temperature correlation
1135
temp_bins = np.arange(-10, 41, 5)
1136
temp_binned = pd.cut(cell_temp, temp_bins)
1137
temp_grouped = combined_losses.groupby(temp_binned)['temp_factor'].mean()
1138

1139
axes[2, 0].plot(temp_grouped.index.map(lambda x: x.mid), temp_grouped, 
1140
               'ro-', linewidth=2, markersize=6)
1141
axes[2, 0].set_xlabel('Cell Temperature (°C)')
1142
axes[2, 0].set_ylabel('Temperature Factor')
1143
axes[2, 0].set_title('Temperature Loss Relationship')
1144
axes[2, 0].grid(True)
1145

1146
# Add theoretical curve
1147
temp_range = np.arange(-5, 40, 1)
1148
theoretical = 1 + temp_coeff * (temp_range - 25)
1149
axes[2, 0].plot(temp_range, theoretical, 'b--', 
1150
               linewidth=2, label='Theoretical')
1151
axes[2, 0].legend()
1152

1153
# AOI correlation 
1154
aoi_bins = np.arange(0, 91, 10)
1155
aoi_binned = pd.cut(aoi, aoi_bins)
1156
aoi_grouped = combined_losses.groupby(aoi_binned)['iam_factor'].mean()
1157

1158
axes[2, 1].plot(aoi_grouped.index.map(lambda x: x.mid), aoi_grouped,
1159
               'go-', linewidth=2, markersize=6)
1160
axes[2, 1].set_xlabel('Angle of Incidence (degrees)')
1161
axes[2, 1].set_ylabel('IAM Factor')
1162
axes[2, 1].set_title('Incidence Angle Loss Relationship')
1163
axes[2, 1].grid(True)
1164

1165
# Monthly energy comparison
1166
monthly_energy = pd.DataFrame({
1167
    'baseline': combined_losses['baseline_poa'].resample('M').sum() / 1000,
1168
    'with_losses': combined_losses['effective_irradiance'].resample('M').sum() / 1000
1169
})
1170
monthly_energy['loss_pct'] = (1 - monthly_energy['with_losses'] / monthly_energy['baseline']) * 100
1171

1172
x_months = range(1, 13)
1173
width = 0.35
1174

1175
bars1 = axes[3, 0].bar([x - width/2 for x in x_months], monthly_energy['baseline'], 
1176
                      width, label='Baseline', alpha=0.7)
1177
bars2 = axes[3, 0].bar([x + width/2 for x in x_months], monthly_energy['with_losses'], 
1178
                      width, label='With Losses', alpha=0.7)
1179

1180
axes[3, 0].set_xlabel('Month')
1181
axes[3, 0].set_ylabel('Monthly Energy (kWh/m²)')
1182
axes[3, 0].set_title('Monthly Energy: Baseline vs With Losses')
1183
axes[3, 0].set_xticks(x_months)
1184
axes[3, 0].legend()
1185
axes[3, 0].grid(True, axis='y')
1186

1187
# System-level annual summary
1188
system_baseline_kwh = baseline_energy * system_params['system_power'] / 1000
1189
system_actual_kwh = effective_energy * system_params['system_power'] / 1000
1190
system_loss_kwh = system_baseline_kwh - system_actual_kwh
1191

1192
economic_impact = system_loss_kwh * 0.12  # Assume $0.12/kWh
1193

1194
# Summary pie chart
1195
loss_breakdown = {
1196
    'Temperature': loss_impacts['temp_factor'],
1197
    'PVWatts Losses': loss_impacts['pvwatts_losses'], 
1198
    'Soiling': loss_impacts['soiling_factor'],
1199
    'Spectral': loss_impacts['spectral_factor'],
1200
    'IAM': loss_impacts['iam_factor'],
1201
    'DC Ohmic': loss_impacts['dc_ohmic_factor'],
1202
    'Degradation': loss_impacts['degradation']
1203
}
1204

1205
# Only show significant losses (>0.5%)
1206
significant_losses = {k: v for k, v in loss_breakdown.items() if v > 0.5}
1207
other_losses = sum(v for k, v in loss_breakdown.items() if v <= 0.5)
1208
if other_losses > 0:
1209
    significant_losses['Other'] = other_losses
1210

1211
axes[3, 1].pie(significant_losses.values(), labels=significant_losses.keys(), 
1212
               autopct='%1.1f%%', startangle=90)
1213
axes[3, 1].set_title('Annual Loss Breakdown')
1214

1215
plt.tight_layout()
1216
plt.show()
1217

1218
# Print comprehensive summary
1219
print(f"\nSystem Performance Summary:")
1220
print("="*40)
1221
print(f"System size: {system_params['system_power']/1000:.1f} kW")
1222
print(f"Baseline annual energy: {system_baseline_kwh:.0f} kWh")
1223
print(f"Actual annual energy: {system_actual_kwh:.0f} kWh")
1224
print(f"Annual energy loss: {system_loss_kwh:.0f} kWh ({total_loss_pct:.1f}%)")
1225
print(f"Economic impact: ${economic_impact:.0f}/year")
1226

1227
print(f"\nPerformance Ratio: {system_actual_kwh/system_baseline_kwh:.3f}")
1228
print(f"Capacity Factor: {system_actual_kwh/(system_params['system_power']*8760/1000):.1f}%")
1229

1230
print(f"\nTop Loss Contributors:")
1231
top_losses = sorted(loss_impacts.items(), key=lambda x: x[1], reverse=True)[:5]
1232
for i, (factor, impact) in enumerate(top_losses, 1):
1233
    print(f"{i}. {factor.replace('_', ' ').title()}: {impact:.2f}%")
1234
```