Full project code

This post contains the full source code for all three MATLAB scripts developed for the PDE3823 final year project: IoT-Based Wind Energy Monitoring and Analysis System.

Matlab File 1: wind_turbine_analysis.m

%% Wind Turbine Data Analysis - Main Script

% Author: Muhammad Cezanne Niyaz

% Date: February 2026

% Description: Analyzes simulated wind turbine dataset and generates plots

clear all;

close all;

clc;

%% 1. Load the Dataset

fprintf('Loading dataset...\n');

data = readtable('sample_motor_dataset.csv');

% Extract columns (handles both with and without underscores)

if ismember('Time_s', data.Properties.VariableNames)

time = data.Time_s;

voltage = data.Voltage_V;

current = data.Current_A;

rpm = data.RPM;

power = data.Power_W;

else

time = data.Time;

voltage = data.Voltage;

current = data.Current;

rpm = data.RPM;

power = data.Power;

end

fprintf('Dataset loaded successfully!\n');

fprintf('Total data points: %d\n', length(time));

fprintf('Time range: %.1f to %.1f seconds\n', min(time), max(time));

fprintf('Voltage range: %.2f to %.2f V\n', min(voltage), max(voltage));

fprintf('Current range: %.3f to %.3f A\n', min(current), max(current));

fprintf('RPM range: %.1f to %.1f RPM\n', min(rpm), max(rpm));

fprintf('Power range: %.3f to %.3f W\n\n', min(power), max(power));

%% 2. Plot 1: Voltage, Current, and Power vs Time

figure(1);

set(gcf, 'Position', [100, 100, 1200, 800]);

% Voltage subplot

subplot(3,1,1);

plot(time, voltage, 'b-', 'LineWidth', 2);

grid on;

xlabel('Time (s)', 'FontSize', 12, 'FontWeight', 'bold');

ylabel('Voltage (V)', 'FontSize', 12, 'FontWeight', 'bold');

title('Voltage vs Time', 'FontSize', 14, 'FontWeight', 'bold');

ylim([0 12]);

% Current subplot

subplot(3,1,2);

plot(time, current, 'r-', 'LineWidth', 2);

grid on;

xlabel('Time (s)', 'FontSize', 12, 'FontWeight', 'bold');

ylabel('Current (A)', 'FontSize', 12, 'FontWeight', 'bold');

title('Current vs Time', 'FontSize', 14, 'FontWeight', 'bold');

ylim([0 0.8]);

% Power subplot

subplot(3,1,3);

plot(time, power, 'g-', 'LineWidth', 2);

grid on;

xlabel('Time (s)', 'FontSize', 12, 'FontWeight', 'bold');

ylabel('Power (W)', 'FontSize', 12, 'FontWeight', 'bold');

title('Power vs Time', 'FontSize', 14, 'FontWeight', 'bold');

ylim([0 5]);

sgtitle('Wind Turbine Performance Over Time', 'FontSize', 16, 'FontWeight', 'bold');

% Save figure

saveas(gcf, 'plot1_time_series.png');

fprintf('Plot 1 saved: plot1_time_series.png\n');

%% 3. Plot 2: RPM vs Power

figure(2);

set(gcf, 'Position', [100, 100, 800, 600]);

scatter(rpm, power, 50, 'filled', 'MarkerFaceColor', [0.2 0.6 0.8]);

hold on;

% Add trend line

p = polyfit(rpm, power, 2); % 2nd order polynomial fit

rpm_fit = linspace(min(rpm), max(rpm), 100);

power_fit = polyval(p, rpm_fit);

plot(rpm_fit, power_fit, 'r-', 'LineWidth', 2);

grid on;

xlabel('RPM', 'FontSize', 12, 'FontWeight', 'bold');

ylabel('Power (W)', 'FontSize', 12, 'FontWeight', 'bold');

title('Power vs RPM', 'FontSize', 14, 'FontWeight', 'bold');

legend('Data Points', 'Polynomial Fit', 'Location', 'northwest');

% Add correlation coefficient

R = corrcoef(rpm, power);

text(min(rpm)+50, max(power)-0.3, sprintf('R = %.3f', R(1,2)), ...

'FontSize', 11, 'BackgroundColor', 'white');

saveas(gcf, 'plot2_rpm_vs_power.png');

fprintf('Plot 2 saved: plot2_rpm_vs_power.png\n');

%% 4. Plot 3: Voltage vs Current

figure(3);

set(gcf, 'Position', [100, 100, 800, 600]);

scatter(voltage, current, 50, power, 'filled');

colorbar;

colormap('jet');

grid on;

xlabel('Voltage (V)', 'FontSize', 12, 'FontWeight', 'bold');

ylabel('Current (A)', 'FontSize', 12, 'FontWeight', 'bold');

title('Current vs Voltage (colored by Power)', 'FontSize', 14, 'FontWeight', 'bold');

clim([0 5]);

saveas(gcf, 'plot3_voltage_vs_current.png');

fprintf('Plot 3 saved: plot3_voltage_vs_current.png\n');

%% 5. Plot 4: Efficiency Analysis

figure(4);

set(gcf, 'Position', [100, 100, 800, 600]);

% Calculate theoretical maximum power (V * I_max)

power_max_theoretical = voltage * 0.8; % Max current is 0.8A

efficiency = (power ./ power_max_theoretical) * 100;

plot(time, efficiency, 'b-', 'LineWidth', 2);

grid on;

xlabel('Time (s)', 'FontSize', 12, 'FontWeight', 'bold');

ylabel('Efficiency (%)', 'FontSize', 12, 'FontWeight', 'bold');

title('System Efficiency Over Time', 'FontSize', 14, 'FontWeight', 'bold');

ylim([0 100]);

% Add mean efficiency line

mean_eff = mean(efficiency);

hold on;

yline(mean_eff, 'r--', 'LineWidth', 2);

legend('Efficiency', sprintf('Mean = %.1f%%', mean_eff), 'Location', 'best');

saveas(gcf, 'plot4_efficiency.png');

fprintf('Plot 4 saved: plot4_efficiency.png\n');

%% 6. Statistical Summary

fprintf('\n=== STATISTICAL SUMMARY ===\n\n');

fprintf('Voltage Statistics:\n');

fprintf(' Mean: %.2f V\n', mean(voltage));

fprintf(' Std Dev: %.2f V\n', std(voltage));

fprintf(' Min: %.2f V, Max: %.2f V\n\n', min(voltage), max(voltage));

fprintf('Current Statistics:\n');

fprintf(' Mean: %.3f A\n', mean(current));

fprintf(' Std Dev: %.3f A\n', std(current));

fprintf(' Min: %.3f A, Max: %.3f A\n\n', min(current), max(current));

fprintf('RPM Statistics:\n');

fprintf(' Mean: %.1f RPM\n', mean(rpm));

fprintf(' Std Dev: %.1f RPM\n', std(rpm));

fprintf(' Min: %.1f RPM, Max: %.1f RPM\n\n', min(rpm), max(rpm));

fprintf('Power Statistics:\n');

fprintf(' Mean: %.3f W\n', mean(power));

fprintf(' Std Dev: %.3f W\n', std(power));

fprintf(' Min: %.3f W, Max: %.3f W\n', min(power), max(power));

fprintf(' Total Energy: %.3f J (over 60s)\n\n', sum(power) * 1); % 1 second intervals

fprintf('Efficiency Statistics:\n');

fprintf(' Mean: %.2f%%\n', mean(efficiency));

fprintf(' Std Dev: %.2f%%\n', std(efficiency));

fprintf(' Min: %.2f%%, Max: %.2f%%\n\n', min(efficiency), max(efficiency));

%% 7. Save processed data

output_table = table(time, voltage, current, rpm, power, efficiency, ...

'VariableNames', {'Time_s', 'Voltage_V', 'Current_A', 'RPM', 'Power_W', 'Efficiency_pct'});

writetable(output_table, 'processed_data.csv');

fprintf('Processed data saved: processed_data.csv\n');

fprintf('\nAnalysis complete!\n');

MATLAB File 2: multi_wind_turbine_second_draft_.m

%% Multi-Turbine Wind Farm Simulation

% Author: Muhammad Cezanne Niyaz

% Date: February 2026

% Description: Simulates multiple turbines with different characteristics

% feeding into a single grid

clear all;

close all;

clc;

%% 1. Define Turbine Parameters

fprintf('=== MULTI-TURBINE WIND FARM SIMULATION ===\n\n');

% Number of turbines

num_turbines = 4;

% Turbine specifications (each turbine has different characteristics)

turbine_specs = struct();

% Turbine 1: Small turbine, low wind

turbine_specs(1).name = 'Turbine 1 (Small, Low Wind)';

turbine_specs(1).diameter = 0.5; % meters

turbine_specs(1).tsr = 3; % Tip Speed Ratio

turbine_specs(1).wind_speed = 3; % m/s

turbine_specs(1).color = [0.2 0.4 0.8];

% Turbine 2: Medium turbine, medium wind

turbine_specs(2).name = 'Turbine 2 (Medium, Med Wind)';

turbine_specs(2).diameter = 0.7;

turbine_specs(2).tsr = 4;

turbine_specs(2).wind_speed = 4.5;

turbine_specs(2).color = [0.8 0.4 0.2];

% Turbine 3: Large turbine, high wind

turbine_specs(3).name = 'Turbine 3 (Large, High Wind)';

turbine_specs(3).diameter = 0.9;

turbine_specs(3).tsr = 5;

turbine_specs(3).wind_speed = 6;

turbine_specs(3).color = [0.2 0.8 0.4];

% Turbine 4: Medium turbine, variable wind

turbine_specs(4).name = 'Turbine 4 (Medium, Variable)';

turbine_specs(4).diameter = 0.6;

turbine_specs(4).tsr = 3.5;

turbine_specs(4).wind_speed = 5;

turbine_specs(4).color = [0.8 0.2 0.8];

%% 2. Calculate RPM for Each Turbine

% Formula: RPM = (60 �� V �� TSR) / (�� �� D)

fprintf('Calculating RPM for each turbine...\n');

for i = 1:num_turbines

V = turbine_specs(i).wind_speed;

TSR = turbine_specs(i).tsr;

D = turbine_specs(i).diameter;

turbine_specs(i).rpm = (60 * V * TSR) / (pi * D);

fprintf('%s: RPM = %.1f\n', turbine_specs(i).name, turbine_specs(i).rpm);

end

%% 3. Simulate Power Output Over Time

time = 0:60; % 0 to 60 seconds

num_samples = length(time);

% Initialize power arrays

power_turbines = zeros(num_turbines, num_samples);

fprintf('\nSimulating power output...\n');

for i = 1:num_turbines

% Base values scaled by turbine characteristics

base_voltage = 4 + (turbine_specs(i).wind_speed / 6) * 8; % 4-12V range

base_power = (turbine_specs(i).wind_speed / 6) * 4; % Scale power with wind speed

for t = 1:num_samples

% Add variation to simulate real wind fluctuations

voltage_var = base_voltage + randn() * 0.8;

voltage_var = max(2, min(12, voltage_var)); % Clamp to 2-12V

% Power varies with wind conditions (sinusoidal + noise)

power_factor = 0.8 + 0.4 * sin(t/10) + randn() * 0.15;

power_factor = max(0.3, min(1.2, power_factor));

% Calculate power for this turbine

power_turbines(i, t) = base_power * power_factor;

power_turbines(i, t) = max(0.5, min(5, power_turbines(i, t))); % Clamp to 0.5-5W

end

end

%% 4. Calculate Total Grid Power

total_power = sum(power_turbines, 1);

%% 5. Plot Individual Turbine Outputs

figure(1);

set(gcf, 'Position', [100, 100, 1200, 800]);

set(gcf, 'Color', 'k'); % Set figure background to black

for i = 1:num_turbines

subplot(2, 2, i);

plot(time, power_turbines(i, :), 'Color', turbine_specs(i).color, 'LineWidth', 2);

% Set black background and white text

set(gca, 'Color', 'k'); % Black background

set(gca, 'XColor', 'w'); % White X-axis

set(gca, 'YColor', 'w'); % White Y-axis

set(gca, 'GridColor', 'w'); % White grid lines

set(gca, 'GridAlpha', 0.3); % Make grid slightly transparent

grid on;

xlabel('Time (s)', 'FontSize', 11, 'FontWeight', 'bold', 'Color', 'w');

ylabel('Power (W)', 'FontSize', 11, 'FontWeight', 'bold', 'Color', 'w');

title(turbine_specs(i).name, 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

xlim([0 60]);

ylim([0 6]);

% Add statistics

mean_power = mean(power_turbines(i, :));

text(2, 5.5, sprintf('Mean: %.2f W\nRPM: %.0f', mean_power, turbine_specs(i).rpm), ...

'FontSize', 9, 'BackgroundColor', 'black', 'Color', 'white', 'EdgeColor', 'white');

end

% Set main title to white

sgtitle('Individual Turbine Power Output', 'FontSize', 14, 'FontWeight', 'bold', 'Color', 'w');

saveas(gcf, 'multi_turbine_individual.png');

fprintf('Plot saved: multi_turbine_individual.png\n');

%% 6. Plot Total Farm Output

figure(2);

set(gcf, 'Position', [100, 100, 1200, 600]);

set(gcf, 'Color', 'k'); % Set figure background to black

% Plot all turbines and total

subplot(1, 2, 1);

hold on;

for i = 1:num_turbines

plot(time, power_turbines(i, :), 'Color', turbine_specs(i).color, ...

'LineWidth', 1.5, 'DisplayName', sprintf('Turbine %d', i));

end

plot(time, total_power, 'w-', 'LineWidth', 3, 'DisplayName', 'Total Grid Power');

% Set black background and white text

set(gca, 'Color', 'k');

set(gca, 'XColor', 'w');

set(gca, 'YColor', 'w');

set(gca, 'GridColor', 'w');

set(gca, 'GridAlpha', 0.3);

grid on;

xlabel('Time (s)', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

ylabel('Power (W)', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

title('All Turbines + Total Grid Power', 'FontSize', 13, 'FontWeight', 'bold', 'Color', 'w');

xlim([0 60]);

legend('Location', 'best', 'TextColor', 'w', 'Color', 'k', 'EdgeColor', 'w');

% Bar chart of average power per turbine

subplot(1, 2, 2);

mean_powers = mean(power_turbines, 2);

bar_colors = zeros(num_turbines, 3);

for i = 1:num_turbines

bar_colors(i, :) = turbine_specs(i).color;

end

bar(1:num_turbines, mean_powers, 'FaceColor', 'flat', 'CData', bar_colors);

% Set black background and white text

set(gca, 'Color', 'k');

set(gca, 'XColor', 'w');

set(gca, 'YColor', 'w');

set(gca, 'GridColor', 'w');

set(gca, 'GridAlpha', 0.3);

grid on;

xlabel('Turbine Number', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

ylabel('Average Power (W)', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

title('Average Power per Turbine', 'FontSize', 13, 'FontWeight', 'bold', 'Color', 'w');

xticks(1:num_turbines);

% Add total power text

hold on;

yline(mean(total_power), 'r--', 'LineWidth', 2);

text(1.5, mean(total_power) + 0.3, sprintf('Grid Mean: %.2f W', mean(total_power)), ...

'FontSize', 11, 'FontWeight', 'bold', 'Color', 'r');

sgtitle('Multi-Turbine Wind Farm Analysis', 'FontSize', 15, 'FontWeight', 'bold', 'Color', 'w');

saveas(gcf, 'multi_turbine_total.png');

fprintf('Plot saved: multi_turbine_total.png\n');

%% 7. Power Distribution Analysis

figure(3);

set(gcf, 'Position', [100, 100, 800, 600]);

set(gcf, 'Color', 'k'); % Set figure background to black

% Create stacked area chart

area(time, power_turbines');

colormap(bar_colors);

% Set black background and white text

set(gca, 'Color', 'k');

set(gca, 'XColor', 'w');

set(gca, 'YColor', 'w');

set(gca, 'GridColor', 'w');

set(gca, 'GridAlpha', 0.3);

grid on;

xlabel('Time (s)', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

ylabel('Power (W)', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

title('Power Contribution by Turbine (Stacked)', 'FontSize', 14, 'FontWeight', 'bold', 'Color', 'w');

xlim([0 60]);

legend(arrayfun(@(i) sprintf('Turbine %d', i), 1:num_turbines, 'UniformOutput', false), ...

'Location', 'northwest', 'TextColor', 'w', 'Color', 'k', 'EdgeColor', 'w');

saveas(gcf, 'multi_turbine_stacked.png');

fprintf('Plot saved: multi_turbine_stacked.png\n');

%% 8. Statistical Summary

fprintf('\n=== MULTI-TURBINE STATISTICS ===\n\n');

for i = 1:num_turbines

fprintf('%s:\n', turbine_specs(i).name);

fprintf(' RPM: %.1f\n', turbine_specs(i).rpm);

fprintf(' Mean Power: %.3f W\n', mean(power_turbines(i, :)));

fprintf(' Std Dev: %.3f W\n', std(power_turbines(i, :)));

fprintf(' Peak Power: %.3f W\n\n', max(power_turbines(i, :)));

end

fprintf('Total Grid Output:\n');

fprintf(' Mean Power: %.3f W\n', mean(total_power));

fprintf(' Std Dev: %.3f W\n', std(total_power));

fprintf(' Peak Power: %.3f W\n', max(total_power));

fprintf(' Total Energy: %.3f J (over 60s)\n\n', sum(total_power) * 1);

%% 9. Identify Problems

fprintf('=== IDENTIFIED PROBLEMS ===\n\n');

% Problem 1: Power Discrepancy

power_range = max(mean_powers) - min(mean_powers);

fprintf('1. Power Discrepancy:\n');

fprintf(' Range: %.3f W (%.1f%% variation)\n', ...

power_range, (power_range / mean(mean_powers)) * 100);

fprintf(' Solution: Implement power aggregation and balancing\n\n');

% Problem 2: Voltage Instability (simulated)

fprintf('2. Voltage Instability:\n');

fprintf(' Different turbines produce different voltages\n');

fprintf(' Solution: DC-DC converters to standardize voltage\n\n');

% Problem 3: Efficiency

total_theoretical = sum(mean_powers) * 1.2; % Could be 20% better

actual_efficiency = (mean(total_power) / total_theoretical) * 100;

fprintf('3. System Efficiency:\n');

fprintf(' Current: %.1f%%\n', actual_efficiency);

fprintf(' Solution: Optimize power management algorithms\n\n');

%% 10. Save Results

results_table = table();

results_table.Time_s = time';

for i = 1:num_turbines

results_table.(sprintf('Turbine%d_W', i)) = power_turbines(i, :)';

end

results_table.Total_Power_W = total_power';

writetable(results_table, 'multi_turbine_results.csv');

fprintf('Results saved: multi_turbine_results.csv\n');

fprintf('\nMulti-turbine analysis complete!\n');

MATLAB File 3: advanced_analysis_multiturbine.m

%% Advanced Multi-Turbine Wind Farm Analysis

% Author: Muhammad Cezanne Niyaz

% Date: March 2026

% Description: Advanced analysis comparing different turbine configurations,

% grid stability issues, and optimization strategies

clear all;

close all;

clc;

fprintf('=== ADVANCED MULTI-TURBINE ANALYSIS ===\n\n');

%% SCENARIO 1: Different Turbine Sizes (Same Wind Speed)

fprintf('SCENARIO 1: Testing different turbine sizes with same wind (5 m/s)\n');

scenario1_turbines = [

0.4, 3.5, 5; % Small: 0.4m diameter, TSR=3.5, 5 m/s wind

0.6, 3.5, 5; % Medium-Small

0.8, 3.5, 5; % Medium-Large

1.0, 3.5, 5; % Large

];

[power1, rpm1, stats1] = simulate_scenario(scenario1_turbines, 'Same Wind, Different Sizes');

%% SCENARIO 2: Different Wind Speeds (Same Turbine Size)

fprintf('\nSCENARIO 2: Testing same turbine size (0.7m) with different wind speeds\n');

scenario2_turbines = [

0.7, 4, 3; % Low wind: 3 m/s

0.7, 4, 4.5; % Medium wind: 4.5 m/s

0.7, 4, 6; % High wind: 6 m/s

0.7, 4, 7.5; % Very high wind: 7.5 m/s

];

[power2, rpm2, stats2] = simulate_scenario(scenario2_turbines, 'Same Size, Different Wind');

%% SCENARIO 3: Optimized Configuration (Balanced)

fprintf('\nSCENARIO 3: Optimized balanced configuration\n');

scenario3_turbines = [

0.6, 4, 4.5; % All turbines similar specs

0.6, 4, 4.5;

0.6, 4, 4.5;

0.6, 4, 4.5;

];

[power3, rpm3, stats3] = simulate_scenario(scenario3_turbines, 'Optimized Balanced');

%% SCENARIO 4: Worst Case (Maximum Variation)

fprintf('\nSCENARIO 4: Worst case - maximum variation\n');

scenario4_turbines = [

0.3, 3, 2; % Very small, very low wind

0.5, 3.5, 4; % Small-medium

0.8, 4.5, 6.5; % Large

1.2, 5, 8; % Very large, very high wind

];

[power4, rpm4, stats4] = simulate_scenario(scenario4_turbines, 'Worst Case Variation');

%% COMPARATIVE ANALYSIS

fprintf('\n=== COMPARATIVE ANALYSIS ===\n\n');

% Compile all scenarios

all_scenarios = {

'Scenario 1: Different Sizes', stats1;

'Scenario 2: Different Winds', stats2;

'Scenario 3: Optimized', stats3;

'Scenario 4: Worst Case', stats4;

};

% Plot comparison

figure(5);

set(gcf, 'Position', [100, 100, 1400, 900]);

set(gcf, 'Color', 'k'); % Black background

% Subplot 1: Power Variation Comparison

subplot(2, 2, 1);

variations = [stats1.power_variation, stats2.power_variation, ...

stats3.power_variation, stats4.power_variation];

bar(variations, 'FaceColor', [0.2 0.6 0.8]);

set(gca, 'Color', 'k');

set(gca, 'XColor', 'w');

set(gca, 'YColor', 'w');

set(gca, 'GridColor', 'w');

set(gca, 'GridAlpha', 0.3);

set(gca, 'XTickLabel', {'Different Sizes', 'Different Winds', 'Optimized', 'Worst Case'});

ylabel('Power Variation (%)', 'FontSize', 11, 'FontWeight', 'bold', 'Color', 'w');

title('Power Variation Across Scenarios', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

grid on;

xtickangle(45);

% Subplot 2: Total Grid Power Comparison

subplot(2, 2, 2);

total_powers = [stats1.total_mean_power, stats2.total_mean_power, ...

stats3.total_mean_power, stats4.total_mean_power];

bar(total_powers, 'FaceColor', [0.8 0.4 0.2]);

set(gca, 'Color', 'k');

set(gca, 'XColor', 'w');

set(gca, 'YColor', 'w');

set(gca, 'GridColor', 'w');

set(gca, 'GridAlpha', 0.3);

set(gca, 'XTickLabel', {'Different Sizes', 'Different Winds', 'Optimized', 'Worst Case'});

ylabel('Mean Grid Power (W)', 'FontSize', 11, 'FontWeight', 'bold', 'Color', 'w');

title('Total Grid Power Output', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

grid on;

xtickangle(45);

% Subplot 3: Grid Stability Score

subplot(2, 2, 3);

stability_scores = [stats1.stability_score, stats2.stability_score, ...

stats3.stability_score, stats4.stability_score];

bar(stability_scores, 'FaceColor', [0.2 0.8 0.4]);

set(gca, 'Color', 'k');

set(gca, 'XColor', 'w');

set(gca, 'YColor', 'w');

set(gca, 'GridColor', 'w');

set(gca, 'GridAlpha', 0.3);

set(gca, 'XTickLabel', {'Different Sizes', 'Different Winds', 'Optimized', 'Worst Case'});

ylabel('Stability Score (0-100)', 'FontSize', 11, 'FontWeight', 'bold', 'Color', 'w');

title('Grid Stability Analysis', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

grid on;

xtickangle(45);

ylim([0 100]);

% Subplot 4: Efficiency Score

subplot(2, 2, 4);

efficiency_scores = [stats1.efficiency, stats2.efficiency, ...

stats3.efficiency, stats4.efficiency];

bar(efficiency_scores, 'FaceColor', [0.8 0.2 0.8]);

set(gca, 'Color', 'k');

set(gca, 'XColor', 'w');

set(gca, 'YColor', 'w');

set(gca, 'GridColor', 'w');

set(gca, 'GridAlpha', 0.3);

set(gca, 'XTickLabel', {'Different Sizes', 'Different Winds', 'Optimized', 'Worst Case'});

ylabel('System Efficiency (%)', 'FontSize', 11, 'FontWeight', 'bold', 'Color', 'w');

title('Overall System Efficiency', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

grid on;

xtickangle(45);

sgtitle('Multi-Turbine Scenario Comparison', 'FontSize', 15, 'FontWeight', 'bold', 'Color', 'w');

saveas(gcf, 'scenario_comparison.png');

fprintf('Saved: scenario_comparison.png\n');

%% POWER VS RPM DETAILED ANALYSIS

fprintf('\n=== POWER VS RPM ANALYSIS ===\n\n');

figure(6);

set(gcf, 'Position', [100, 100, 1200, 800]);

set(gcf, 'Color', 'k'); % Black background

colors_plot = [0.2 0.4 0.8; 0.8 0.4 0.2; 0.2 0.8 0.4; 0.8 0.2 0.8];

% Flatten power matrices properly for each scenario

power1_flat = power1(:);

power2_flat = power2(:);

power3_flat = power3(:);

power4_flat = power4(:);

% Create RPM arrays that match the power data size

rpm1_expanded = repmat(rpm1, 1, size(power1, 2));

rpm1_flat = rpm1_expanded(:);

rpm2_expanded = repmat(rpm2, 1, size(power2, 2));

rpm2_flat = rpm2_expanded(:);

rpm3_expanded = repmat(rpm3, 1, size(power3, 2));

rpm3_flat = rpm3_expanded(:);

rpm4_expanded = repmat(rpm4, 1, size(power4, 2));

rpm4_flat = rpm4_expanded(:);

% Combine all scenario data

all_rpm = [rpm1_flat; rpm2_flat; rpm3_flat; rpm4_flat];

all_power = [power1_flat; power2_flat; power3_flat; power4_flat];

% Plot 1: Combined view

subplot(2, 2, 1);

scatter(all_rpm, all_power, 20, 'filled', 'MarkerFaceAlpha', 0.6);

hold on;

% Polynomial fit

p = polyfit(all_rpm, all_power, 2);

rpm_fit = linspace(min(all_rpm), max(all_rpm), 100);

power_fit = polyval(p, rpm_fit);

plot(rpm_fit, power_fit, 'r-', 'LineWidth', 2);

set(gca, 'Color', 'k');

set(gca, 'XColor', 'w');

set(gca, 'YColor', 'w');

set(gca, 'GridColor', 'w');

set(gca, 'GridAlpha', 0.3);

xlabel('RPM', 'FontSize', 11, 'FontWeight', 'bold', 'Color', 'w');

ylabel('Power (W)', 'FontSize', 11, 'FontWeight', 'bold', 'Color', 'w');

title('Combined Power vs RPM (All Scenarios)', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

grid on;

legend('Data Points', 'Polynomial Fit', 'Location', 'northwest', 'TextColor', 'w', 'Color', 'k', 'EdgeColor', 'w');

% Plot 2: Scenario 1

subplot(2, 2, 2);

scatter(rpm1_flat, power1_flat, 30, colors_plot(1,:), 'filled');

set(gca, 'Color', 'k');

set(gca, 'XColor', 'w');

set(gca, 'YColor', 'w');

set(gca, 'GridColor', 'w');

set(gca, 'GridAlpha', 0.3);

xlabel('RPM', 'FontSize', 11, 'FontWeight', 'bold', 'Color', 'w');

ylabel('Power (W)', 'FontSize', 11, 'FontWeight', 'bold', 'Color', 'w');

title('Scenario 1: Different Sizes', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

grid on;

% Plot 3: Scenario 2

subplot(2, 2, 3);

scatter(rpm2_flat, power2_flat, 30, colors_plot(2,:), 'filled');

set(gca, 'Color', 'k');

set(gca, 'XColor', 'w');

set(gca, 'YColor', 'w');

set(gca, 'GridColor', 'w');

set(gca, 'GridAlpha', 0.3);

xlabel('RPM', 'FontSize', 11, 'FontWeight', 'bold', 'Color', 'w');

ylabel('Power (W)', 'FontSize', 11, 'FontWeight', 'bold', 'Color', 'w');

title('Scenario 2: Different Winds', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

grid on;

% Plot 4: Scenario 3

subplot(2, 2, 4);

scatter(rpm3_flat, power3_flat, 30, colors_plot(3,:), 'filled');

set(gca, 'Color', 'k');

set(gca, 'XColor', 'w');

set(gca, 'YColor', 'w');

set(gca, 'GridColor', 'w');

set(gca, 'GridAlpha', 0.3);

xlabel('RPM', 'FontSize', 11, 'FontWeight', 'bold', 'Color', 'w');

ylabel('Power (W)', 'FontSize', 11, 'FontWeight', 'bold', 'Color', 'w');

title('Scenario 3: Optimized', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

grid on;

sgtitle('Power vs RPM Analysis by Scenario', 'FontSize', 15, 'FontWeight', 'bold', 'Color', 'w');

saveas(gcf, 'power_vs_rpm_detailed.png');

fprintf('Saved: power_vs_rpm_detailed.png\n');

%% GRID STABILITY HEATMAP

fprintf('\n=== GRID STABILITY HEATMAP ===\n\n');

figure(7);

set(gcf, 'Position', [100, 100, 1000, 700]);

set(gcf, 'Color', 'k'); % Black background

% Create stability matrix

stability_matrix = [

stats1.power_variation, stats1.voltage_std, stats1.total_mean_power;

stats2.power_variation, stats2.voltage_std, stats2.total_mean_power;

stats3.power_variation, stats3.voltage_std, stats3.total_mean_power;

stats4.power_variation, stats4.voltage_std, stats4.total_mean_power;

];

% Normalize for heatmap

stability_matrix_norm = stability_matrix ./ max(stability_matrix);

imagesc(stability_matrix_norm);

colormap(flipud(hot));

colorbar;

set(gca, 'Color', 'k');

set(gca, 'XColor', 'w');

set(gca, 'YColor', 'w');

set(gca, 'XTick', 1:3);

set(gca, 'XTickLabel', {'Power Variation', 'Voltage Std Dev', 'Total Power'});

set(gca, 'YTick', 1:4);

set(gca, 'YTickLabel', {'Different Sizes', 'Different Winds', 'Optimized', 'Worst Case'});

xtickangle(45);

title('Grid Stability Metrics (Normalized)', 'FontSize', 14, 'FontWeight', 'bold', 'Color', 'w');

xlabel('Metric', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

ylabel('Scenario', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'w');

% Add text values

for i = 1:4

for j = 1:3

text(j, i, sprintf('%.2f', stability_matrix_norm(i,j)), ...

'HorizontalAlignment', 'center', 'Color', 'white', ...

'FontWeight', 'bold', 'FontSize', 11);

end

end

saveas(gcf, 'stability_heatmap.png');

fprintf('Saved: stability_heatmap.png\n');

%% RECOMMENDATIONS SUMMARY

fprintf('\n=== ANALYSIS RECOMMENDATIONS ===\n\n');

fprintf('Best Configuration: %s\n', all_scenarios{3,1});

fprintf(' - Lowest power variation: %.1f%%\n', stats3.power_variation);

fprintf(' - Highest stability score: %.1f\n', stats3.stability_score);

fprintf(' - Good efficiency: %.1f%%\n\n', stats3.efficiency);

fprintf('Worst Configuration: %s\n', all_scenarios{4,1});

fprintf(' - Highest power variation: %.1f%%\n', stats4.power_variation);

fprintf(' - Lowest stability score: %.1f\n', stats4.stability_score);

fprintf(' - Lower efficiency: %.1f%%\n\n', stats4.efficiency);

fprintf('Key Findings:\n');

fprintf('1. Balanced turbine sizes improve grid stability by %.1f%%\n', ...

((stats3.stability_score - stats1.stability_score) / stats1.stability_score) * 100);

fprintf('2. Power variation reduced from %.1f%% to %.1f%% with optimization\n', ...

stats4.power_variation, stats3.power_variation);

fprintf('3. Total power output maximized in Scenario 2: %.2f W\n', stats2.total_mean_power);

fprintf('\n=== ADVANCED ANALYSIS COMPLETE ===\n');

%% HELPER FUNCTION: Simulate Scenario

function [power_matrix, rpm_array, stats] = simulate_scenario(turbine_config, scenario_name)

% turbine_config: [diameter, tsr, wind_speed] for each turbine

num_turbines = size(turbine_config, 1);

time = 0:60;

num_samples = length(time);

power_matrix = zeros(num_turbines, num_samples);

rpm_array = zeros(num_turbines, 1);

voltage_matrix = zeros(num_turbines, num_samples);

for i = 1:num_turbines

D = turbine_config(i, 1);

TSR = turbine_config(i, 2);

V = turbine_config(i, 3);

% Calculate RPM

rpm_array(i) = (60 * V * TSR) / (pi * D);

% Calculate base power

base_power = (V / 8) * 4; % Scale with wind speed

% Simulate power over time

for t = 1:num_samples

power_factor = 0.8 + 0.4 * sin(t/10) + randn() * 0.15;

power_factor = max(0.3, min(1.2, power_factor));

power_matrix(i, t) = base_power * power_factor;

power_matrix(i, t) = max(0.3, min(6, power_matrix(i, t)));

% Simulate voltage

voltage_matrix(i, t) = 4 + (V / 8) * 8 + randn() * 0.5;

voltage_matrix(i, t) = max(2, min(12, voltage_matrix(i, t)));

end

end

% Calculate statistics

mean_powers = mean(power_matrix, 2);

total_power = sum(power_matrix, 1);

stats.power_variation = (max(mean_powers) - min(mean_powers)) / mean(mean_powers) * 100;

stats.total_mean_power = mean(total_power);

stats.total_std_power = std(total_power);

stats.voltage_std = mean(std(voltage_matrix, 0, 2));

% Stability score (0-100, higher is better)

% Lower variation = higher stability

stats.stability_score = 100 - min(stats.power_variation, 100);

% Efficiency (theoretical max vs actual)

theoretical_max = sum(mean_powers) * 1.2;

stats.efficiency = (stats.total_mean_power / theoretical_max) * 100;

fprintf('%s:\n', scenario_name);

fprintf(' RPM range: %.0f - %.0f\n', min(rpm_array), max(rpm_array));

fprintf(' Power variation: %.1f%%\n', stats.power_variation);

fprintf(' Total mean power: %.2f W\n', stats.total_mean_power);

fprintf(' Stability score: %.1f\n', stats.stability_score);

fprintf(' Efficiency: %.1f%%\n', stats.efficiency);

end