It has been quite some time since I've done surface fitting. Are you looking for an illustration like the one shown below? It is easier for me to visualize the relationship between two inputs and one output. I think it is possible for you to group the clusters according to the number of houses.
Update: I think that the general equation for the grid energy would probably look like a nested polynomial function:
where 
, 
, 
are third-order polynomial functions given by
, , are third-order polynomial functions given by
.%% Data for 10 Houses
figure(1)
x10 = [12e3 12e3 12e3 24e3 24e3 24e3 32e3 32e3 32e3]'; % PV Capacity
y10 = [ 3e4 6e4 10e4 3e4 6e4 10e4 3e4 6e4 10e4]'; % Battery
z10 = [37240 23430 5840 15980 698 698.4 12810 0 0]'; % Grid Energy
plot3(x10, y10, z10, 'rp'), grid on,
hold on
% Surface (manual fitting)
XX = linspace(x10(1), x10(end), 51);
YY = linspace(y10(1), y10(end), 51);
[X, Y] = meshgrid(XX, YY);
% z(x, y) = f(x)*y^2 + g(x)*y + h(x)
f10 = (-3.645e-14)*X.^2 + 1.894e-09*X - 1.719e-05;
g10 = (3.997e-09)*X.^2 - 0.0002003*X + 1.341;
h10 = (-1.842e-05)*X.^2 + 0.06133*X + 5.35e+04;
Z10 = f10.*Y.^2 + g10.*Y + h10;
s = surf(X, Y, Z10, 'FaceAlpha', 0.25); s.EdgeColor = 'none';
view(145, 30)
hold off
% Labels
xlabel('PV Capacity / (Watt)'), ylabel('Battery / (Wh)'), zlabel('Grid Energy / (Wh)')
title('Number of Houses: 10')
%% Data for 20 Houses
figure(2)
x20 = [24e3 24e3 24e3 48e3 48e3 48e3 64e3 64e3 64e3]'; % PV Capacity
y20 = [ 3e4 6e4 10e4 3e4 6e4 10e4 3e4 6e4 10e4]'; % Battery
z20 = [93680 78670 57030 61840 35960 8150 56370 32280 7308]'; % Grid Energy
plot3(x20, y20, z20, 'rp'), grid on,
hold on
% Surface (manual fitting)
XX = linspace(x20(1), x20(end), 51);
YY = linspace(y20(1), y20(end), 51);
[X, Y] = meshgrid(XX, YY);
% z(x, y) = f(x)*y^2 + g(x)*y + h(x)
f20 = (-2.845e-15)*X.^2 + (3.287e-10)*X - 6.832e-06;
g20 = (7.266e-10)*X.^2 + (-7.856e-05)*X + 1.019;
h20 = (5.354e-06)*X.^2 + (-1.035)*X + 1.293e+05;
Z20 = f20.*Y.^2 + g20.*Y + h20;
s = surf(X, Y, Z20, 'FaceAlpha', 0.25); s.EdgeColor = 'none';
view(145, 30)
hold off
% Labels
xlabel('PV Capacity / (Watt)'), ylabel('Battery / (Wh)'), zlabel('Grid Energy / (Wh)')
title('Number of Houses: 20')
%% Data for 30 Houses
figure(3)
x30 = [36e3 36e3 36e3 72e3 72e3 72e3 96e3 96e3 96e3]'; % PV Capacity
y30 = [ 3e4 6e4 10e4 3e4 6e4 10e4 3e4 6e4 10e4]'; % Battery
z30 = [120500 106600 90520 106400 74190 46690 93280 75360 42610]'; % Grid Energy
plot3(x30, y30, z30, 'rp'), grid on,
hold on
% Surface (manual fitting)
XX = linspace(x30(1), x30(end), 51);
YY = linspace(y30(1), y30(end), 51);
[X, Y] = meshgrid(XX, YY);
% z(x, y) = f(x)*y^2 + g(x)*y + h(x)
f30 = (-8.176e-15)*X.^2 + 1.012e-09*X - 2.496e-05;
g30 = (1.349e-09)*X.^2 - 0.0001742*X + 3.982;
h30 = (-3.565e-05)*X.^2 + 4.197*X + 3.11e+04;
Z30 = f30.*Y.^2 + g30.*Y + h30;
s = surf(X, Y, Z30, 'FaceAlpha', 0.25); s.EdgeColor = 'none';
view(145, 30)
hold off
% Labels
xlabel('PV Capacity / (Watt)'), ylabel('Battery / (Wh)'), zlabel('Grid Energy / (Wh)')
title('Number of Houses: 30')
%% Data for 50 Houses
x50 = [6e4 6e4 6e4 12e4 12e4 12e4 16e4 16e4 16e4]'; % PV Capacity
y50 = [3e4 6e4 10e4 3e4 6e4 10e4 3e4 6e4 10e4]'; % Battery
z50 = [279300 262100 211600 193100 186400 131100 187200 164200 133500]'; % Grid Energy
plot3(x50, y50, z50, 'rp'), grid on,
hold on
% Surface (manual fitting)
XX = linspace(x50(1), x50(end), 51);
YY = linspace(y50(1), y50(end), 51);
[X, Y] = meshgrid(XX, YY);
% z(x, y) = f(x)*y^2 + g(x)*y + h(x)
f50 = (5.256e-15)*X.^2 - 1.058e-09*X + 3.471e-05;
g50 = (-6.672e-10)*X.^2 + 0.000136*X - 5.445;
h50 = (2.818e-05)*X.^2 - 6.886*X + 5.905e+05;
Z50 = f50.*Y.^2 + g50.*Y + h50;
s = surf(X, Y, Z50, 'FaceAlpha', 0.25); s.EdgeColor = 'none';
view(145, 30)
hold off
% Labels
xlabel('PV Capacity / (Watt)'), ylabel('Battery / (Wh)'), zlabel('Grid Energy / (Wh)')
title('Number of Houses: 50')
