How to plot a multivariable function with constraints

Hello,

I plotted a 2-variables function and I would like to highlight (or plot/intersect) only the z points for which x (k in the function) and y satisfy the following constraints: 250020*((pi/4)*k^2*y)-325.33=0

I have been trying several approaches but I cannot make it happen. I hope you can help.

Thank you!

f=@(k,y) (8120*(((pi*((((k+(2*(5*0.0026))+(2*0.0026))+(((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350))*2)+((4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))*2)).^2)-((((k+(2*(5*0.0026))+(2*0.0026))+(((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350))*2)+((4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))*2))-(2*((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350)))).^2))*y)/4)+(((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))*(4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))*y.*24)))+(8120*(pi*((k.^2)-((k-(2*((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350)))).^2))*y)/4)+(7450*((pi*y.*(((k+(2*(5*0.0026))).^2)-(k.^2)))/4)+(8933*pi*(0.003/2).^2*((15*(2*4))*((2*y)+((2*(pi*((k+(2*(5*0.0026))+(2*0.0026))+(4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350)))))))/8)+(4*(((pi*((k+(2*(5*0.0026))+(2*0.0026))+(4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))))*3));
fsurf(f);
xlim([0.02 0.1]);
ylim([0.04 0.4]);
zlim([0 100]);

Answer 1

[John D'Errico]

Can you solve for y, as a function of k, given the constraint? (YES.)

If so, then do you really have a TWO variable problem? (NO.)

syms k y
ysol = solve(250020*((pi/4)*k^2*y)-325.33==0,y)

ysol = 

z = (8120*(((pi*((((k+(2*(5*0.0026))+(2*0.0026))+(((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350))*2)+((4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))*2)).^2)-((((k+(2*(5*0.0026))+(2*0.0026))+(((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350))*2)+((4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))*2))-(2*((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350)))).^2))*y)/4)+(((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))*(4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))*y.*24)))+(8120*(pi*((k.^2)-((k-(2*((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350)))).^2))*y)/4)+(7450*((pi*y.*(((k+(2*(5*0.0026))).^2)-(k.^2)))/4)+(8933*pi*(0.003/2).^2*((15*(2*4))*((2*y)+((2*(pi*((k+(2*(5*0.0026))+(2*0.0026))+(4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350)))))))/8)+(4*(((pi*((k+(2*(5*0.0026))+(2*0.0026))+(4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))))*3));
z_k = subs(z,y,ysol)

z_k = 

It may look messy, but who cares?

fplot(z_k,[0.02,0.1])
xlabel 'k'
ylabel 'z'

Again, this is not a 2 variable problem, due to the constraint.

Can you plot it in 3-dimensions anyway? Well, yes.

K = linspace(0.02,0.1,100);
yfun = matlabFunction(ysol);
Y = yfun(K);
f = @(k,y) (8120*(((pi*((((k+(2*(5*0.0026))+(2*0.0026))+(((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350))*2)+((4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))*2)).^2)-((((k+(2*(5*0.0026))+(2*0.0026))+(((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350))*2)+((4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))*2))-(2*((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350)))).^2))*y)/4)+(((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))*(4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))*y.*24)))+(8120*(pi*((k.^2)-((k-(2*((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350)))).^2))*y)/4)+(7450*((pi*y.*(((k+(2*(5*0.0026))).^2)-(k.^2)))/4)+(8933*pi*(0.003/2).^2*((15*(2*4))*((2*y)+((2*(pi*((k+(2*(5*0.0026))+(2*0.0026))+(4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350)))))))/8)+(4*(((pi*((k+(2*(5*0.0026))+(2*0.0026))+(4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))))*3));
f = str2func(vectorize(func2str(f)));
plot3(K,Y,f(K,Y),'o-')
grid on
box on
xlabel 'k'
ylabel 'y'
zlabel 'z'

[Lorenzo Manni]

Thank you very much!

This does the job perfectly. However, I realised that the very long function does not behave as expected, probably due to a typo.

How could I use the "short version" below?

y = (rhoL*(((pi*((Dso^2)-((Dso-(2*lsy))^2))*y)/4)+(tw*Sd*y*Nsl)))+(rhoL*(pi*((k^2)-((k-(2*lsy))^2))*y)/4)+(rhoM*((pi*y*((Dmo^2)-(k^2)))/4)+(rhoC*pi*(Dstra/2)^2*lc*Nph));

Where everything is defined as:

%% Motor Requirements

Trwheel = 3090.625;                     %           torque at rear wheels
Nwheel = 2152.820;                      %           max wheel speed
wwheel = 225.443;                       %           max wheel speed

GR = 9.5;                               %           gear ratio (9.5 or 7.56)

Vdc = 1000;                             % [V]       max pattery DC voltage
Pout = 250000;                          % [W]       max power out of Battery
Tmotor = Trwheel/GR;                    % [Nm]      motor torque required
Nsmax = Nwheel*GR;                      % [rpm]     max synchonous speed
wmax = wwheel*GR;                       % [rad/s]   max synchonous speed
Nctp = (60/(2*pi))*Pout/Trwheel;        % [rpm]     Wheel Speed at which Power and Torque are Max (transition from constant T to contsant P region)
Nsbase = Nctp*GR;                       % [rpm]     base speed
wbase = ((2*pi)/60)*Nsbase;             % [rpm]     base speed
opT = 60;                               % [°C]      operating temperature
Nph = 3;                                %           number of phases
%% Material Properties

% Wire (Copper)
T0 = 20;                                % [°C]      reference temperature
rhoCu0 = 1.7241e-8;                     % [Ωm]      Copper resistivity
alpha = 0.004;                          %           thermal resisitvity coefficient
rhoCures = rhoCu0*(1+(alpha*(opT-T0))); % [Ωm]      Copper resistivity at operating T
Dstra = 3e-4;                           % [m]       diameter of each conductor strand 
rhoC = 8933;                            % [kg/m^3]  Copper density

% Lamination (VX48, 49% Co-Fe)
Bs = 2.350;                             % [T]       saturation flux density
tauL = 1e-4;                            % [m]       lamination thickness
rhores = 4.2e-7;                        % [Ωm]      lamination resistivity
sigmaL = 1/rhores;                      % [S/m]     lamination conductivity
rhoL = 8120;                            % [kg/m^3]  lamination density
PFE = 1.5;                              % [W/kg]    lamination loss density @50Hz, 1.5T
Stckf = 0.97;                           %           stacking factor [Chau, 2015 p85]

% Magnets (N55, Nd-Fe-B)
Br = 1.49;                              % [T]       remeance magnetic flux density
rhoM = 7450;                            % [kg/m^3]  magnet density

%% Motor Architecture [INTPUT]

% Sizing Parameters 
p = 8;                                  %           poles
Nsl = 24;                               %           slots 
gcs = 0.5e-3;                           % [m]       gap between carbon sleeve and stator (Zhao, 2019)
% Sd = 48.47e-3;                        % [m]       slot depth
Nc = 15;                                %           turns per coil
SFF = 0.6;                              %           slot fill factor (MSO 0.9 CHin, 2021)
JCu = 30e6;                             % [A/m^2]   RMS max current density (Chu, 2020)
pf = 1;                                 %           percentage gap between poles

% Carbon Sleeve
tcs = 0.002;        %LOOP?              % [m]       carbon sleeve thickness (Golonov, 2018)
lg = tcs + gcs;                         % [m]       air gap length between stator and rotor
YScs = 800e6;                           % [N/m^2]   yield strength of carbon fibre
SF = 1.4;                               %           safety factor
rhoCS = 1750;                           % [kg/m^3]  carbon sleeve density (Blisset, 2019, p145)

%% Motor Evaluation [COMPUTATION]

% General Sizing
pp = p/2;                               %           pole pairs               
lm = 5*lg;                              % [m]       magnet thickness
f = (pp*Nsbase)/60;                     % [Hz]      line voltage frequency

% Diameters as a function of k
Dmo = k+(2*lm);                       % [m]       magnets outer diameter 
Dsb = k+(2*lm)+(2*lg);                % [m]       stator bore diameter ESTIMATE

% Air Gap Flux
Bg = (Br*lm)/(lm+lg);                   % [T]       air gap magnetic flux density
Am = (pi*Dsb*y)/(p/pf);  %CHECK L       % [m^2]     single pole area at stator bore

% Stator (Tooth, Yoke, Slot Width & Outer Diameter)
phim = Am*Bg;                           % [Tm^2]    flux per pole per unit length
phit = phim/(Nsl);   %CHECK             % [Tm^2]    flux per tooth
tw = (2*pi()*(Dmo/2)*Bg)/(Nsl*Stckf*Bs);% [m]       tooth width
tpitch = (pi*Dsb)/Nsl;                  % [m]       stator tooth pitch
Sw = tpitch-tw;                         % [m]       slot width (0.45 ~ 0.62)*tpitch [Chau, 2015 p63]
% Sw = tpitch*0.535;
Sd = 4.5*Sw;                            % [m]       slot depth (3.5 ~ 5.5)*bss [Chau, 2015 p63]
lsy = (pi*(Dmo/2)*Bg)/(p*Stckf*Bs);     % [m]       stator & rotor yoke thickness
Ds = Dsb+Sd;                            % [m]       Davg in middle of slot
Swa = ((pi*Ds)/Nsl)-tw;                 % [m]       average width of the slot (trapezoidal)
Dso = Dsb+(lsy*2)+(Sd*2);               % [m]       stator outer diameter
Yr = (lsy)/(lsy+Sd);                    %           yoke ratio

% Current
As = Swa*Sd;                            % [m^2]     slot C-S area
ACu = SFF*As;                           % [m^2]     copper C-S area per slot
Irms = ACu*JCu;                         % [A]       RMS current in each slot

% Conductor Current
slp = Nsl/(p*Nph);                     %           slots per pole per phase
Z = 2*Nc;                               %           total number of conductors per slot
i = Irms/Z;                             % [A]       peak line current

% Wire Area
Acon = ACu/Z;                           % [m^2]     conductor C-S area per slot
Dcon = sqrt((4*Acon)/pi);               % [m]       conductor diameter
Astra = pi*((Dstra/2)^2);               % [m^2]     area of each strand
Nstra = Acon/Astra;                     %           number of strands

% MLT
Wc = (pi*Ds)/p;                         % [m]       coil width
MLT = (2*y)+((2*(pi*Ds))/p)+(4*Swa);    % [m]       mean length of one turn

% Resistance
Ncoil = 2*pp;                           %           number of coils
Ntph = Nc*Ncoil;                        %           total no. of turns in series per phase
lc = Ntph*MLT;                          % [m]       total length of wire per phase
R = (rhoCures*lc)/Acon;                 % [Ohm]     resistance per phase
PCu = Nph*(i^2)*R;                      % [W]       copper losses

The problem I encounter is that if I use the function with all the named terms, when I run the f = @(k,y) obviously it doesn't recognise them.

Many thanks

[John D'Errico]

I think you need to think about what are variables. Can you just put the name k into an equation without ever having defined it before? For example, I see these lines:

% Diameters as a function of k
Dmo = k+(2*lm);                       % [m]       magnets outer diameter 

Unrecognized function or variable 'k'.

Dsb = k+(2*lm)+(2*lg);                % [m]       stator bore diameter ESTIMATE

Was k ever defined? NO. How should matlab treat those lines? It will say this:

Unrecognized function or variable 'k'. 

And how about y? What is y? Yet you use it without any predefinition...

% Air Gap Flux
...
Am = (pi*Dsb*y)/(p/pf);  %CHECK L       % [m^2]     single pole area at stator bore

So perhaps you need to look more carefully at what I did. Did I use syms there? What does syms do? READ THE HELP! Learn MATLAB, rather than just typing your equations and hoping MATLAB will understand how to interpret what you write. You may sometimes get lucky, but not often.

What I did was to define k and y as symbolic unknowns.

help syms

doc syms

If you do that, then all of your variables now become parameterized by the unknowns k and y. So now you might try this (although I have no clue what you re really doing in this, what is your goal):

fimplicit(y == (rhoL*(((pi*((Dso^2)-((Dso-(2*lsy))^2))*y)/4)+(tw*Sd*y*Nsl)))+(rhoL*(pi*((k^2)-((k-(2*lsy))^2))*y)/4)+(rhoM*((pi*y*((Dmo^2)-(k^2)))/4)+(rhoC*pi*(Dstra/2)^2*lc*Nph)))

That is because the expression you write is an implicit function of k and y.

[Lorenzo Manni]

Thank you. I will try to be more clear.

I've built an analytical parametrised model of a SPM electric motor. All the equations use values from the material property section and the sizing parameters section. From these, I developed an equation to calculate the mass of the motor.

Then, I've set the variables k (stator diameter) and y (active length), so I could ran in a separate code fmincon to evaluate the minimum mass achievable with a combination of k and y that would satisfy the constraints discussed in the main post (which refer to the torque output as a function of the diameter k and length y).

I needed to plot the overall effect of k and y on the mass (the surface) and then highlight the combinations that satisfy the torque requirement (what you helped me with).

The issue is that, to generate the mass function I posted initially, I did manually input all the different parameters until the only terms remaining were the variables k and y. However, I would like to run the function with different combination of parameters without having to change them manually (i.e. assign a different rhoL).

This is what I have so far:

%% Motor Requirements

Trwheel = 3090.625;                     %           torque at rear wheels
Nwheel = 2152.820;                      %           max wheel speed
wwheel = 225.443;                       %           max wheel speed

GR = 9.5;                               %           gear ratio (9.5 or 7.56)

Vdc = 1000;                             % [V]       max pattery DC voltage
Pout = 250000;                          % [W]       max power out of Battery
Tmotor = Trwheel/GR;                    % [Nm]      motor torque required
Nsmax = Nwheel*GR;                      % [rpm]     max synchonous speed
wmax = wwheel*GR;                       % [rad/s]   max synchonous speed
Nctp = (60/(2*pi))*Pout/Trwheel;        % [rpm]     Wheel Speed at which Power and Torque are Max (transition from constant T to contsant P region)
Nsbase = Nctp*GR;                       % [rpm]     base speed
wbase = ((2*pi)/60)*Nsbase;             % [rpm]     base speed
opT = 60;                               % [°C]      operating temperature
Nph = 3;                                %           number of phases

%% Material Properties

% Wire (Copper)
T0 = 20;                                % [°C]      reference temperature
rhoCu0 = 1.7241e-8;                     % [Ωm]      Copper resistivity
alpha = 0.004;                          %           thermal resisitvity coefficient
rhoCures = rhoCu0*(1+(alpha*(opT-T0))); % [Ωm]      Copper resistivity at operating T
Dstra = 3e-4;                           % [m]       diameter of each conductor strand 
rhoC = 8933;                            % [kg/m^3]  Copper density

% Lamination (VX48, 49% Co-Fe)
Bs = 2.350;                             % [T]       saturation flux density
tauL = 1e-4;                            % [m]       lamination thickness
rhores = 4.2e-7;                        % [Ωm]      lamination resistivity
sigmaL = 1/rhores;                      % [S/m]     lamination conductivity
rhoL = 8120;                            % [kg/m^3]  lamination density
PFE = 1.5;                              % [W/kg]    lamination loss density @50Hz, 1.5T
Stckf = 0.97;                           %           stacking factor [Chau, 2015 p85]

% Magnets (N55, Nd-Fe-B)
Br = 1.49;                              % [T]       remeance magnetic flux density
rhoM = 7450;                            % [kg/m^3]  magnet density

%% Motor Architecture [INTPUT]

% Sizing Parameters 
p = 8;                                  %           poles
Nsl = 24;                               %           slots 
gcs = 0.5e-3;                           % [m]       gap between carbon sleeve and stator (Zhao, 2019)
% Sd = 48.47e-3;                        % [m]       slot depth
Nc = 15;                                %           turns per coil
SFF = 0.6;                              %           slot fill factor (MSO 0.9 CHin, 2021)
JCu = 30e6;                             % [A/m^2]   RMS max current density (Chu, 2020)
pf = 1;                                 %           percentage gap between poles

% Carbon Sleeve
tcs = 0.002;        %LOOP?              % [m]       carbon sleeve thickness (Golonov, 2018)
lg = tcs + gcs;                         % [m]       air gap length between stator and rotor
YScs = 800e6;                           % [N/m^2]   yield strength of carbon fibre
SF = 1.4;                               %           safety factor
rhoCS = 1750;                           % [kg/m^3]  carbon sleeve density (Blisset, 2019, p145)

%% Motor Evaluation [COMPUTATION]

% General Sizing
pp = p/2;                               %           pole pairs               
lm = 5*lg;                              % [m]       magnet thickness
f = (pp*Nsbase)/60;                     % [Hz]      line voltage frequency

% Diameters as a function of k
Dmo = k+(2*lm);                       % [m]       magnets outer diameter 
Dsb = k+(2*lm)+(2*lg);                % [m]       stator bore diameter ESTIMATE

% Air Gap Flux
Bg = (Br*lm)/(lm+lg);                   % [T]       air gap magnetic flux density
Am = (pi*Dsb*y)/(p/pf);  %CHECK L       % [m^2]     single pole area at stator bore

% Stator (Tooth, Yoke, Slot Width & Outer Diameter)
phim = Am*Bg;                           % [Tm^2]    flux per pole per unit length
phit = phim/(Nsl);   %CHECK             % [Tm^2]    flux per tooth
tw = (2*pi()*(Dmo/2)*Bg)/(Nsl*Stckf*Bs);% [m]       tooth width
tpitch = (pi*Dsb)/Nsl;                  % [m]       stator tooth pitch
Sw = tpitch-tw;                         % [m]       slot width (0.45 ~ 0.62)*tpitch [Chau, 2015 p63]
% Sw = tpitch*0.535;
Sd = 4.5*Sw;                            % [m]       slot depth (3.5 ~ 5.5)*bss [Chau, 2015 p63]
lsy = (pi*(Dmo/2)*Bg)/(p*Stckf*Bs);     % [m]       stator & rotor yoke thickness
Ds = Dsb+Sd;                            % [m]       Davg in middle of slot
Swa = ((pi*Ds)/Nsl)-tw;                 % [m]       average width of the slot (trapezoidal)
Dso = Dsb+(lsy*2)+(Sd*2);               % [m]       stator outer diameter
Yr = (lsy)/(lsy+Sd);                    %           yoke ratio

% Current
As = Swa*Sd;                            % [m^2]     slot C-S area
ACu = SFF*As;                           % [m^2]     copper C-S area per slot
Irms = ACu*JCu;                         % [A]       RMS current in each slot

% Conductor Current
slp = Nsl/(p*Nph);                     %           slots per pole per phase
Z = 2*Nc;                               %           total number of conductors per slot
i = Irms/Z;                             % [A]       peak line current

% Wire Area
Acon = ACu/Z;                           % [m^2]     conductor C-S area per slot
Dcon = sqrt((4*Acon)/pi);               % [m]       conductor diameter
Astra = pi*((Dstra/2)^2);               % [m^2]     area of each strand
Nstra = Acon/Astra;                     %           number of strands

% MLT
Wc = (pi*Ds)/p;                         % [m]       coil width
MLT = (2*y)+((2*(pi*Ds))/p)+(4*Swa);    % [m]       mean length of one turn

% Resistance
Ncoil = 2*pp;                           %           number of coils
Ntph = Nc*Ncoil;                        %           total no. of turns in series per phase
lc = Ntph*MLT;                          % [m]       total length of wire per phase
R = (rhoCures*lc)/Acon;                 % [Ohm]     resistance per phase
PCu = Nph*(i^2)*R;                      % [W]       copper losses

%% Plotting

% Surface Plot

ysol = solve(250020*((pi/4)*k^2*y)-325.33==0,y);
% z = (8120*(((pi*((((k+(2*(5*0.0026))+(2*0.0026))+(((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350))*2)+((4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))*2)).^2)-((((k+(2*(5*0.0026))+(2*0.0026))+(((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350))*2)+((4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))*2))-(2*((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350)))).^2))*y)/4)+(((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))*(4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))*y.*24)))+(8120*(pi*((k.^2)-((k-(2*((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350)))).^2))*y)/4)+(7450*((pi*y.*(((k+(2*(5*0.0026))).^2)-(k.^2)))/4)+(8933*pi*(0.003/2).^2*((15*(2*4))*((2*y)+((2*(pi*((k+(2*(5*0.0026))+(2*0.0026))+(4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350)))))))/8)+(4*(((pi*((k+(2*(5*0.0026))+(2*0.0026))+(4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))))*3));
z = (rhoL*(((pi*((Dso^2)-((Dso-(2*lsy))^2))*y)/4)+(tw*Sd*y*Nsl)))+(rhoL*(pi*((k^2)-((k-(2*lsy))^2))*y)/4)+(rhoM*((pi*y*((Dmo^2)-(k^2)))/4)+(rhoC*pi*(Dstra/2)^2*lc*Nph));
z_k = subs(z,y,ysol);
% fplot(z_k,[0.02,0.1]);
% xlabel 'k';
% ylabel 'z';
K = linspace(0.02,0.1,100);
yfun = matlabFunction(ysol);
Y = yfun(K);
% f = @(k,y) (8120*(((pi*((((k+(2*(5*0.0026))+(2*0.0026))+(((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350))*2)+((4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))*2)).^2)-((((k+(2*(5*0.0026))+(2*0.0026))+(((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350))*2)+((4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))*2))-(2*((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350)))).^2))*y)/4)+(((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))*(4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))*y.*24)))+(8120*(pi*((k.^2)-((k-(2*((pi*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(8*0.97*2.350)))).^2))*y)/4)+(7450*((pi*y.*(((k+(2*(5*0.0026))).^2)-(k.^2)))/4)+(8933*pi*(0.003/2).^2*((15*(2*4))*((2*y)+((2*(pi*((k+(2*(5*0.0026))+(2*0.0026))+(4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350)))))))/8)+(4*(((pi*((k+(2*(5*0.0026))+(2*0.0026))+(4.5*(((pi*(k+(2*(5*0.0026))+(2*0.0026)))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))))/24)-((2*pi()*((k+(2*(5*0.0026)))/2)*((1.49*(5*0.0026))/((5*0.0026)+0.0026)))/(24*0.97*2.350))))))*3));
% fimplicit(z == (rhoL*(((pi*((Dso^2)-((Dso-(2*lsy))^2))*y)/4)+(tw*Sd*y*Nsl)))+(rhoL*(pi*((k^2)-((k-(2*lsy))^2))*y)/4)+(rhoM*((pi*y*((Dmo^2)-(k^2)))/4)+(rhoC*pi*(Dstra/2)^2*lc*Nph)))
fsurf(z)
% grid on
% box on
xlim([0.02 0.1]);
ylim([0.04 0.4]);
zlim([0 80]);
ax = gca;
ax.XAxis.Exponent = -3;
ax.YAxis.Exponent = -3;
xlabel('Stator Diameter (m)');
ylabel('Active Length (m)');
zlabel('Motor Mass (kg)');
hold on
f = str2func(vectorize(func2str(z)));
plot3(K,Y,z(K,Y),'r-','LineWidth',2.0)
hold off

It works until I plot the surface. However I cannot plot the line you did. How could I fix it?

Many thanks for your help

[Lorenzo Manni]

The issue is that now the line

f = str2func(vectorize(func2str(f)))

does not work because there is no

f=@(k,y) 

before, since I am writing the function using

syms k y

How would you fix that @John D'Errico?

Many thanks