Solve a system of equations

Can anyone help me to write a code for solving the following system of equations?

w0=10; k1=3; k2=6; V1=20; V2=30;

x21=(w0+x13)/(w0+x13+k1*V1);

x22=(w0+x13);

x31=x22*x21/(x22+k2*V2);

x32=x22;

x13=x31*x32;

Suppose that x13 has an initial value (any number, let's say it is 1); Then we define what x21, x22, x31, 32 will be equal to and new value of x13*n;

If |x13*n-x|>=e (where e is a small number), then a new value of x13 should be x13*n, and we should solve that system until |x13*n-x13|<e;

Below is my idea, but I know it may be stupid but don't judge too harsh because I'm new to matlab and don't know how to do that.

w0=10; k1=3; k2=6; V1=20; V2=30; x13=1; e=0.00001; n=0.1
while abs(x13*n-x13)>=e
x21=(w0+x13)/(w0+x13+k1*V1);
x22=w0+x13;
x31=x22*x21/(x22+k2*V2);
x32=x22;
x13=x31*x32;
x13=x13*n;
end 

It never stops counting

[darova]

Not allowed to use fsolve()?

[Rustem Devletov]

We can use anything. If you can, show me please how to do that

Answer 1

[Torsten]

x22 - w0 =
x13 = 
x31 * x32 =
x22 * x21 / (x22+k2*V2) * x22 =
x22 * (x22/(x22+k1*V1)) / (x22+k2*V2) * x22

This gives a quadratic equation in x22. Solve it.

Once you have x22, x13 = x22 - w0.

[Rustem Devletov]

I didn't understand... What I have to do is to give an initial guess to x13

[Torsten]

x22 - w0 = x22 * (x22/(x22+k1*V1)) / (x22+k2*V2) * x22

This follows from your relations. Solve this equation for x22. Since x22 = x13 + w0, it follows that x13 = x22 - w0. So you don't need any nonlinear solver or iteration.

[Rustem Devletov]

Hm..alright. I'll try to

[Rustem Devletov]

file 1

function F=solvx22(y)
w0=10; k1=3; k2=6; V1=20; V2=30; x13=1; e=0.001; n=0.1;
F=y - w0 - y * (y/(y+k1*V1)) / (y+k2*V2) * y
end

file 2

fun=@solvx22
x0=1;
x=fsolve(fun,x0)

Is that what you mean?

[Torsten]

Can't you rewrite the equation

x22 - w0 = x22 * (x22/(x22+k1*V1)) / (x22+k2*V2) * x22

as

a*x22^2 + b*x22 + c = 0

by calculation with pencil and paper and solve for x22 ?

[Rustem Devletov]

AAahhhh okaaay, now I understand what you want. But this is the thing, I have to solve it with iterations

[Torsten]

Then try your fsolve-code. fsolve uses Newton's (iterative) method. But it might converge to the wrong solution because the quadratic equation in x22 usually has two solutions.

[Rustem Devletov]

Can you check it out please?

file 1

global  w0 k1 k2 x32 
w0=10; k1=3; k2=6; V1=20; V2=30 x13=1; e=0.00001;
k=x13;
x21=(w0+x13)/(w0+x13+k1*V1);
x22=w0+x13;
x31=x22*x21/(x22+k2*V2);
x32=x22;
x13=x31*x32;

while abs(x13-k)>=e
k=x13;
x21=(w0+x13)/(w0+x13+k1*V1);
x22=w0+x13;
x31=x22*x21/(x22+k2*V2);
x32=x22;
x13=x31*x32
end 

file 2

function [Q]= myfun(x)
global x32 
alpha1=0.5;
alpha2=0.5;
V1=x(1); 
V2=x(2);
C1=10; %cost 1
C2=15; %cost 2
C3=24;%cost 3 
Q=C1*V1^alpha1+C2*V2*alpha2+C3*x32;
end

file 3 (solution)

clear all
clc
fun = @myfun;
x = [20, 30];
A = [];
b = [];
Aeq = [];
beq = [];
lb = [10, 20];
ub = [30, 40];
nonlcon = [];
options = optimoptions('fmincon');
Cany 
[x,fval,exitflag,output] = fmincon(fun,x,A,b,Aeq,beq,lb,ub,nonlcon,options)

I get the following error

Error using fmincon (line 612)
Supplied objective function must return a scalar value.

Error in solution (line 14)
[x,fval,exitflag,output] = fmincon(fun,x,A,b,Aeq,beq,lb,ub,nonlcon,options)