0×1 empty double column vector

My Code -

disp("Break Away Point")
syms sigma omega k f real
s = sigma+j*omega;
fx = (s^3)+(3*s^2)+(2*s)+k;
f = expand(fx);
f1 = real(f);
f2 = imag(f);
A = jacobian([f1,f2],[sigma omega]);
omega = 0;
sigma  = vpasolve(det(eval(A)),sigma)
A2 = eval(f1);
k = double(vpasolve(A2,k))

k =

0×1 empty double column vector

any suggestion how to fix it?

[KSSV]

What you want to do here?

A2 = eval(f1);

[Matt J]

There's nothing obviously invalid about the answer. Maybe the solution set is truly empty.

[Tenzing Thiley]

no it is not if I type A(1) it is giving value

Answer 1

[Walter Roberson]

There are multiple solutions, but your code was trying to find a single k value that satisfied all of the possibilites at the same time.

Also, do not eval() a symbolic expression. eval() of a symbolic expression gives an undefined result.

disp("Break Away Point")

Break Away Point

syms sigma omega k f real
s = sigma+j*omega;
fx = (s^3)+(3*s^2)+(2*s)+k;
f = expand(fx);
f1 = real(f)

f1 = 

f2 = imag(f)

f2 = 

A = jacobian([f1,f2],[sigma omega])

A = 

Omega = 0;
Sigma  = solve(det(subs(A, omega, Omega)),sigma)

Sigma = 

A2 = subs(f1, {sigma, omega}, {Sigma, Omega})

A2 = 

k = arrayfun(@(expr) solve(expr,k), A2)

k = 

[Tenzing Thiley]

How do I set condition to print only value of sigma for which k>0?

[Walter Roberson]

mask = imag(k) == 0 & k > 0;
kpos = k(mask)
sigma_kpos = sigma(mask)

[Tenzing Thiley]

error -The logical indices contain a true value outside of the array bounds.

[Walter Roberson]

mask = imag(k) == 0 & k > 0;
kpos = k(mask)
sigma_kpos = Sigma(mask)

[Tenzing Thiley]

Thank you, I am really great full by quick response. Can you fix one of my problem?

[Tenzing Thiley]

I am trying to convert ss2tf() but getting error. I need SISO output such as u(s)/etta(s), alpha(s)/etta(s) , q(s)/etta(s),theta(s)/etta(s),u(s)/delta(s), alpha(s)/delta(s) , q(s)/delta(s), theta(s)/delta(s).

[Walter Roberson]

Never eval() anything symbolic: eval() is not defined for symbolic expressions or functions.

%% Declaration with parameter
I = 205113.07; % kg-m^2
m= 15118.35;   % kg
c = 3.511;     % m
rho = 1.225;   % kg/m^3
Tm = 49817.6;  % N
S = 37.16;     % m^2
%% Equation of motion(Longitudinal modes)
syms alpha etta delta theta q V 
f1 = ((etta*Tm*cos(alpha))/m) - ((0.5*rho*(V^2)*S*((0.0013*alpha^2)-(0.00438*alpha)+0.1423))./m)-(9.81*sin(theta-alpha));
f2 = q-(etta*Tm*sin(alpha)/(m*V))-(0.5*rho*V*S*((0.0751*alpha)+(0.0144*delta)+0.732)/m)+(9.81*cos(theta-alpha)/V);
f3 = (0.5*rho*V^2*S*c*(-0.00437*alpha-0.0196*delta-0.123*q-0.1885)/I);
f4 = q;
%% Linearization
A = jacobian([f1,f2,f3,f4],[V alpha q theta]);
B = jacobian([f1,f2,f3,f4],[etta delta]);
%% Equillibrium
V = 94.5077;  %m/s
etta = 0.58;
delta = -0.1678; %rad
theta = 0;
alpha = 0;
q = 0;
%% Evaluate System and Control Matrix
format shortG; 
fprintf("State Matrix = ");

State Matrix = 

A = double(subs(A))

A = 4×4

-0.040493 9.8689 0 -9.81 -0.0021967 -0.030908 1 0 -0.013639 -0.015207 -0.42801 0 0 0 1 0

fprintf("Control Matrix = ");

Control Matrix = 

B = double(subs(B))

B = 4×2

3.2952 0 0 -0.0020488 0 -0.068204 0 0

[Tenzing Thiley]

Can you provide me a detail how should I plot a complex graph where I need to use the value of sigma0, poles and sigma1.

clear all;
clc;
%% finding roots
syms sigma omega k f real
s = sigma+j*omega;
ch = 'Characteristics Equation = ';
fx = input(ch);
f = expand(fx);
f1 = real(f);
f2 = imag(f);
k = 0;
omega = 0;
poles = solve(subs(f1) == 0)
if (subs(f2)~=0) 
    A2 = solve(subs(f2) == 0);
end
%% Finding gain for system to be Marginally Stable
syms sigma omega k f real
disp("System is unstable for value of K = ");
sigma = 0;
G = vpasolve(subs(f1),subs(f2),[omega k]);
omega0 = G.omega
k0 = G.k
%% finding Break Away Point
syms sigma omega k f real
disp("Break Away Point");
A = jacobian([f1,f2],[sigma omega]);
Omega = 0;
Sigma  = double(solve(det(subs(A, omega, Omega)),sigma));
A2 = subs(f1, {sigma, omega}, {Sigma, Omega});
k = double(arrayfun(@(expr) solve(expr,k), A2));
mask = imag(k) == 0 & k > 0;
kpos = k(mask);
sigma_kpos = Sigma(mask);
k1 = unique(kpos,'rows')
sigma1 = unique(sigma_kpos,'rows')

[Walter Roberson]

... Don't use clear? Don't overwrite the variables in the first place?

I cannot test your code without knowing the input expected for fx (have you considered writing a function that takes fx as a parameter instead of using input?)

%% finding roots
syms sigma omega k f real
s = sigma+j*omega;
ch = 'Characteristics Equation = ';
fx = input(ch);

Unable to run the 'fevalJSON' function because it calls the 'input' function, which is not supported for this product offering.

f = expand(fx);
f1 = real(f);
f2 = imag(f);
K = 0;
Omega = 0;
f1subs = subs(f1, {k, omega}, {K, Omega})
f2subs = subs(f2, {k, omega}, {K, Omega})
A1 = solve(f1subs == 0)
if (subs(f2)~=0) 
    A2 = solve(f2subs == 0)
end

%% finding unstable gain
disp("System is unstable for value of K = ");
Sigma = 0;
f1subs = subs(f1, {k, omega, sigma}, {K, Omega, Sigma})
f2subs = subs(f2, {k, omega, sigma}, {K, Omega, Sigma})

G = vpasolve([f1subs,f2subs],[omega k])
omega = G.omega
k = G.k
%% finding point

A = jacobian([f1,f2],[sigma omega]);
Omega = 0;
Sigma  = double(solve(det(subs(A, omega, Omega)),sigma));
A2 = subs(f1, {sigma, omega}, {Sigma, Omega})
k = double(arrayfun(@(expr) solve(expr,k), A2))
mask = imag(k) == 0 & k > 0;
kpos = k(mask);
sigma_kpos = Sigma(mask);
k = unique(kpos,'rows')
sigma = unique(sigma_kpos,'rows')

[Tenzing Thiley]

%% finding roots
syms sigma omega k f real
s = sigma+j*omega;
fx = s^3+3*s^2+2*s+k;
f = expand(fx);
f1 = real(f);
f2 = imag(f);
k = 0;
omega = 0;
poles = solve(subs(f1) == 0)
if (subs(f2)~=0) 
    A2 = solve(subs(f2) == 0);
end
%% Finding gain for system to be Marginally Stable
syms sigma omega k f real
disp("System is unstable for value of K = ");
sigma = 0;
G = vpasolve(subs(f1),subs(f2),[omega k]);
omega0 = G.omega
k0 = G.k
%% finding Break Away Point
syms sigma omega k f real
disp("Break Away Point");
A = jacobian([f1,f2],[sigma omega]);
Omega = 0;
Sigma  = double(solve(det(subs(A, omega, Omega)),sigma));
A2 = subs(f1, {sigma, omega}, {Sigma, Omega});
k = double(arrayfun(@(expr) solve(expr,k), A2));
mask = imag(k) == 0 & k > 0;
kpos = k(mask);
sigma_kpos = Sigma(mask);
k1 = unique(kpos,'rows')
sigma1 = unique(sigma_kpos,'rows')

[Tenzing Thiley]

I want to plot a complex graph between sigma and omega and locate point sigma1,omega0,poles but getting error

[Tenzing Thiley]

I treid to use clear K sigma in next code but get error so I overwrite them.

[Walter Roberson]

%% finding roots
syms sigma omega k f real
s = sigma+j*omega;
fx = s^3+3*s^2+2*s+k;
f = expand(fx);
f1 = real(f)

f1 = 

f2 = imag(f)

f2 = 

K = 0;
Omega = 0;
f1subs = subs(f1, {k, omega}, {K, Omega})

f1subs = 

f2subs = subs(f2, {k, omega}, {K, Omega})

f2subs = 

0

poles = solve(f1subs == 0)

poles = 

A2 = sym([])

 
A2 =
 
[ empty sym ]
 

if ~isAlways(f2subs == 0, 'Unknown', 'false')
    A2 = solve(f2subs == 0)
end

%% Finding gain for system to be Marginally Stable
disp("System is unstable for value of K = ");

System is unstable for value of K = 

Sigma = 0;
f1subs = subs(f1, {k, omega, sigma}, {K, Omega, Sigma})

f1subs = 

0

f2subs = subs(f2, {k, omega, sigma}, {K, Omega, Sigma})

f2subs = 

0

G = vpasolve( [f1subs, f2subs], [omega k])

G = struct with fields:
    omega: 0
        k: 0

omega0 = G.omega

omega0 = 

0

k0 = G.k

k0 = 

0

%% finding Break Away Point
disp("Break Away Point");

Break Away Point

A = jacobian([f1,f2],[sigma omega]);
Omega = 0;
eigensigma = solve(det(subs(A, omega, Omega)),sigma)

eigensigma = 

Sigma  = double(eigensigma);
A2 = subs(f1, {sigma, omega}, {Sigma, Omega});
k = double(arrayfun(@(expr) solve(expr,k), A2))

k = 4×1

-0.3849 -0.3849 0.3849 0.3849

mask = imag(k) == 0 & k > 0;
kpos = k(mask);
sigma_kpos = Sigma(mask);
k1 = unique(kpos,'rows')

k1 = 0.3849

sigma1 = unique(sigma_kpos,'rows')

sigma1 = -0.4226

[Tenzing Thiley]

plot?

[Tenzing Thiley]

omega0 and k0 is wrong it should be omega0 =

0

-1.4142135623730950488016887242097

1.4142135623730950488016887242097

k0 =

0

6.0

6.0

[Walter Roberson]

%% finding roots
syms sigma omega k f real
s = sigma+j*omega;
fx = s^3+3*s^2+2*s+k;
f = expand(fx);
f1 = real(f)

f1 = 

f2 = imag(f)

f2 = 

K = 0;
Omega = 0;
f1subs = subs(f1, {k, omega}, {K, Omega})

f1subs = 

f2subs = subs(f2, {k, omega}, {K, Omega})

f2subs = 

0

poles = solve(f1subs == 0)

poles = 

A2 = sym([])

 
A2 =
 
[ empty sym ]
 

if ~isAlways(f2subs == 0, 'Unknown', 'false')
    A2 = solve(f2subs == 0)
end

%% Finding gain for system to be Marginally Stable
disp("System is unstable for value of K = ");

System is unstable for value of K = 

Sigma = 0;
f1subs = subs(f1, {sigma}, {Sigma})

f1subs = 

f2subs = subs(f2, {sigma}, {Sigma})

f2subs = 

G = vpasolve( [f1subs, f2subs], [omega k])

G = struct with fields:
    omega: [3×1 sym]
        k: [3×1 sym]

omega0 = G.omega

omega0 = 

k0 = G.k

k0 = 

%% finding Break Away Point
disp("Break Away Point");

Break Away Point

A = jacobian([f1,f2],[sigma omega]);
Omega = 0;
eigensigma = solve(det(subs(A, omega, Omega)),sigma)

eigensigma = 

Sigma  = double(eigensigma);
A2 = subs(f1, {sigma, omega}, {Sigma, Omega});
k = double(arrayfun(@(expr) solve(expr,k), A2))

k = 4×1

-0.3849 -0.3849 0.3849 0.3849

mask = imag(k) == 0 & k > 0;
kpos = k(mask);
sigma_kpos = Sigma(mask);
k1 = unique(kpos,'rows')

k1 = 0.3849

sigma1 = unique(sigma_kpos,'rows')

sigma1 = -0.4226

[Tenzing Thiley]

plot?

[Walter Roberson]

I want to plot a complex graph between sigma and omega and locate point sigma1,omega0,poles but getting error

You have three variables: sigma, omega, and k. Which variable or variables are to be your independent variable(s) ? Which is to be your dependent variable?

You define fx as an equation, but I do not see any definition of omega in terms of sigma (or the other way around) ? I do not presently see how you would solve for omega in terms of sigma, or for sigma in terms of omega, in order to make one of them the dependent variable. As far as I can see at the moment, sigma and omega are both independent variables.

[Tenzing Thiley]

Graph in complex plane where reult value are

poles

omega0

sigma1

[Tenzing Thiley]

only show this points in complex plane where sigma in x and omega in y

[Walter Roberson]

Potential independent variables:

• real part of sigma
• imaginary part of sigma
• real part of omega
• imaginary part of omega
• k

That seems to be up to 5 independent variables.

[Tenzing Thiley]

Like these value we already obtain from the code only need to locate in complex plane

[Tenzing Thiley]

more clearer pic

[Walter Roberson]

You can use scatter(), and you can use text() if you want to label them.

I notice now that you intend your sigma and omega to be the components of the complex variable s .

But I do not see how k fits into this plot ?

[Tenzing Thiley]

Yeah I ues it earlier but got some error can you show me how you use scatter to obtain plot

[Tenzing Thiley]

Also if answer in complex conjugate it is difficlut to distinguish and we do not need to fit K in plane

[Walter Roberson]

x = [poles; sigma1; zeros(numel(omega0),1)];
y = [zeros(numel(poles),1); 0; omega0(:)];
scatter(x, y, 'x')

[Tenzing Thiley]

Thank You

[Tenzing Thiley]

%% finding roots
syms sigma omega k f real
s = sigma+1i*omega;
fx = s^2*(s+3.6)+(k*(s+0.4));
f = expand(fx);
f = expand(fx);
f1 = real(f);
f2 = imag(f);
K = 0;
Omega = 0;
f1subs = subs(f1, {k, omega}, {K, Omega});
f2subs = subs(f2, {k, omega}, {K, Omega});
K = 0;
Omega = 0;
f1subs = subs(f1, {k, omega}, {K, Omega});
poles = solve(f1subs == 0);
double(poles)
A2 = sym([]);
if ~isAlways(f2subs == 0, 'Unknown', 'false')
    A2 = solve(f2subs == 0);
end
%% Finding gain for system to be Marginally Stable
disp("System is unstable for value of K = ");
Sigma = 0;
f1subs = subs(f1, {sigma}, {Sigma});
f2subs = subs(f2, {sigma}, {Sigma});
G = vpasolve( [f1subs, f2subs], [omega k]);
k0 = unique(G.k,'rows');
disp(k0);
disp("System is unstable for value of omega = ");
omega0 = G.omega;
disp(omega0);
%% finding Break Away Point
A = jacobian([f1,f2],[sigma omega]);
Omega = 0;
eigensigma = solve(det(subs(A, omega, Omega)),sigma);
Sigma  = double(eigensigma);
A2 = subs(f1, {sigma, omega}, {Sigma, Omega});
k = double(arrayfun(@(expr) solve(expr,k), A2));
mask = imag(k) == 0 & k > 0;
kpos = k(mask);
sigma_kpos = Sigma(mask);
k1 = unique(kpos,'rows');
disp("Break Away Point");
sigma1 = unique(sigma_kpos,'rows');
disp(sigma1);

[Tenzing Thiley]

Error using InputOutputModel/feedback (line 137)

The first and second arguments of the "feedback"

command must have compatible I/O sizes.

A = [-0.040493       9.8689            0        -9.81;
    -0.0021967    -0.030908            1            0;
    -0.013639    -0.015207     -0.42801        0;
    0            0            1            0];
B =      [3.2952            0;
    0   -0.0020488;
    0    -0.068204;
    0            0];
C = [1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1];
D = [0 0;0 0;0 0;0 0];
hzpk = zpk(ss(A,B,C,D))
pole(hzpk);
t = [0:0.01:10];
step(0.2*hzpk,t);
axis([0 10 0 0.8]);
ylabel('pitch angle (rad)');
title('Open-loop Step Response');
sys_cl = feedback(hzpk,1);
step(0.2*sys_cl);
ylabel('pitch angle (rad)');
title('Closed-loop Step Response');

[Tenzing Thiley]

Can you check this I am getting error?

[Tenzing Thiley]

?

[Walter Roberson]

You are trying to use the "sign" parameter of feedback(). The "sign" parameter can be added after all of the other positional parameters. However, all of the forms of feedback() require that you pass in two different systems.

[Tenzing Thiley]

here instead of step() I want to input ramp() but the function is not working any idea how I use it?

[Tenzing Thiley]

In above I want to use ramp() input instead of step( ) but the code is giving error as unrecognised function.

can you tell me how should I proceed ?

[Walter Roberson]

slope = randn()

slope = -0.6418

ramp = @(t) (t >= 0) .* t .* slope

ramp = function_handle with value:
    @(t)(t>=0).*t.*slope

fplot(ramp, [-2 10])

In some cases you want the ramp to start from a different location,

starttime = rand() * 5

starttime = 1.8358

ramp2 = @(t) (t >= starttime) .* (t - starttime) .* slope

ramp2 = function_handle with value:
    @(t)(t>=starttime).*(t-starttime).*slope

fplot(ramp2, [-2 10])

[Tenzing Thiley]

But here I have Input hzpk transfer functon

earlier it was easy as step function only require transfer function and time

can you tell me how shoud I proceed with transferfunction as ramp input

[Walter Roberson]

https://www.mathworks.com/matlabcentral/answers/63552-how-find-ramp-response#answer_75146

[Tenzing Thiley]

Do you know how should I use tiledlayout for multiple graph in one page as in this cases I am getting error

Error using DynamicSystem/rlocus (line 65)

This functionality is not supported in TiledChartLayout

%root locus

tiledlayout(2,2)

nexttile

pitch1 = -pitch/(s+10);

rlocus(pitch1);

v = [-4 4 -4 4]; axis(v)

grid

title('Root-Locus Plot of Pitch displacement Auto pilot]');

xlabel('Real Axis')

ylabel('Imag Axis')

%Bode plot

k = 45;

pitch2 = k*pitch1;

margin(pitch2)

grid on;

title('Bode Plot Pitch of displacement Auto pilot');

tiledlayout(1,2)

sys_cl = feedback(pitch2,1);

step(0.2*sys_cl), grid

ylabel('pitch angle (rad)');

title('Closed-loop Step Response')

[Walter Roberson]

You could try using subplot() instead of tiled layout. However, some of the plots in that toolbox need a complete figure to draw in, because the plots they draw are interactive plots instead of static plots