Lösen eines nicht-linearen Gleichungssystems

WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

Günnsen

Forum-Anfänger


	Beiträge: 21

	Anmeldedatum: 26.06.16

	Wohnort: ---

	Version: ---

Verfasst am: 01.08.2016, 11:13 Titel: Lösen eines nicht-linearen Gleichungssystems

Hallo Comunity,

ich arbeite gerade an der Berechnung von Preisen in Abhängigkeit eines Kundenwahlmodells (MNL-Modell). Die Kaufwahrscheinlichkeit für ein Produkt p1 ist abhängig von den Preisen der anderen angeboten Produkte (p2, p3, usw.).

Code:

clear all
syms p1 p2 p3 positive
theta = rand(1); %customer sensitivity to quality
my = rand(1); % random variable
quality = [100, 50, 25]; %quality of each product ordered non increasing
prices = [p1, p2, p3]; % prices

%this formula is used to calculate the buyprobability of a product
buy_prob = (exp((theta.*quality-prices)/my)/(1+ sum(exp(theta.*quality-prices/my))))';

deriv_mat= jacobian(buy_prob); % Jacobi-Matrix
x0=[0;0;0]; %start parameter

Inverse_Matrix = inv(deriv_mat); %Inv of jacobi matrix
h = Inverse_Matrix * (- buy_prob) ;
eqn_sys = gen_eqns (h,p1,p2,p3); %generate equations
[ x,fval ] = solving(x0);

Funktion ohne Link?

h ist ein Vektor der sich wie auf dem Screenshot berechnet. Wobei alpha(p) die Formel für die Kaufwahrscheinlichkeit ist.

Der optimale Preis berechnet sich wie auf dem 2. Screenshot. Da ich den 2. Summanden null setze(also zunächst nicht betrachte), erhalte ich p_t(x)=h(p_t(x)).
Also ein Gleichungssystem.

Mit folgenden Funktionen versuche ich es zu lösen, komme aber nicht weiter.

Code:

function F = gen_eqns( h, p1, p2, p3 )
%GEN_EQNS Summary of this function goes here
% Detailed explanation goes here

F=[h(1)-p1;
h(2)-p2;
h(3)-p3];
end

Funktion ohne Link?

Code:

function [ x,fval ] = solving(x0 )

options = optimset('Display','iter');
[x,fval]= fsolve(@gen_eqns,x0,options);

end

Funktion ohne Link?

Ich habe es mittels einiger Beispiele, dich ich im Netz gefunden habe so versucht.

Komme aber nicht mehr weiter.

Ich hoffe mir kann da jemand weiterhelfen.

Vielen Dank!
Günther

Bildschirmfoto 2016-08-01 um 11.03.32.png

Beschreibung:

Download

Dateiname:

Bildschirmfoto 2016-08-01 um 11.03.32.png

Dateigröße:

51.85 KB

Heruntergeladen:

450 mal

Harald

Forum-Meister


	Beiträge: 24.501

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 01.08.2016, 11:21 Titel:

Hallo,

du hast ja vermutlich Fehlermeldungen bekommen. Poste diese bitte auch.

Du mischt symbolisch und numerisch. Entweder symbolisch (syms + solve) oder numerisch (fsolve).

Angesichts der Komplexität der Gleichungen dürfte eine symbolische Lösung ausscheiden, also numerisch. Schau dir dazu bitte die Hilfe von fsolve und die Beispiele darin an.

Grüße,
Harald

Günnsen

Themenstarter

Forum-Anfänger


	Beiträge: 21

	Anmeldedatum: 26.06.16

	Wohnort: ---

	Version: ---

Verfasst am: 01.08.2016, 11:35 Titel:

Hallo Harald,

hier die Fehlermeldung:

Code:

Error using gen_eqns (line 5)
Not enough input arguments.

Error in fsolve (line 219)
fuser = feval(funfcn{3},x,varargin{:});

Error in solving (line 4)
[x,fval]= fsolve(@gen_eqns,x0,options);

Error in test (line 18)
[ x,fval ] = solving(x0);

Caused by:
Failure in initial user-supplied objective function evaluation. FSOLVE cannot continue.

Funktion ohne Link?

Wenn ich das numerisch löse, darf ich ja keine syms verwenden nehme ich an. Muss ich den Variablen p1,p2,p3 dann einen Wertebereich zuweisen? z.b. 1:0.1:1000 ?

Grüße,
Günther

Harald

Forum-Meister


	Beiträge: 24.501

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 01.08.2016, 14:09 Titel:

Hallo,

nein, du bräuchtest eine Funktion

Code:

function F = equations(p)
p1 = p(1);
p2 = p(2);
p3 = p(3);
% Hier werden aus p1-p3 alle Gleichungen erzeugt

Funktion ohne Link?

Grüße,
Harald

Günnsen

Themenstarter

Forum-Anfänger


	Beiträge: 21

	Anmeldedatum: 26.06.16

	Wohnort: ---

	Version: ---

Verfasst am: 03.08.2016, 16:14 Titel:

Hallo Harald,

vielen Dank für deinen Input. Ich habe es weiterhin symbolisch versucht,
da es mir nicht gelungen ist, eine Jacobimatrix numerisch zu bestimmen.
Nachdem ich mir weitere Tutorials angeschaut habe, bin ich zu folgendem Code gelangt.

Code:

theta = rand(1); %customer sensitivity to quality

my = rand(1);
quality = [100; 50; 25]; %quality of each product ordered non increasing
prices = [x1;x2;x3];
buy_prob =(exp((theta*quality-prices)/my)/(1+ sum(exp(theta*quality-prices/my))));
jacobi = jacobian(buy_prob);
Inverse_Matrix = inv(jacobi);
eqs_sys = prices == Inverse_Matrix * (-1 * buy_prob);
optimal = solve([eqs_sys],[x1, x2, x3]);
x1value = double (optimal.x1);
x2value = double (optimal.y2);
x3value= double (optimal.z3);

Funktion ohne Link?

Der erzeugt aber leider einen Fehler:

Code:

Error using mupadengine/feval (line 163)
The row index is out of range.

Error in solve (line 349)
varargout{1}.(char(vars(i))) = transpose(eng.feval('map', solutions, '_index', i));

Error in test (line 17)
blah = solve([eqs_sys],[x1, x2, x3]);

Funktion ohne Link?

Habe ich etwas nicht bedacht, oder is das Gleichungssystem nicht lösbar?

Grüße,
Günther

Günnsen

Themenstarter

Forum-Anfänger


	Beiträge: 21

	Anmeldedatum: 26.06.16

	Wohnort: ---

	Version: ---

Verfasst am: 03.08.2016, 17:01 Titel:

Hallo,

mit dem Aufruf der Funktion vpasolve erhalte ich nun eine Lösung.

Komisch ist nur , dass ich für x1, x2 und x3 den gleichen Wert erhalte.

*edit
Bei der Berechnung der rechten Seite meines GS sind die Spaltenvektoren identisch.

Harald

Forum-Meister


	Beiträge: 24.501

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 03.08.2016, 19:10 Titel:

Hallo,

setze doch mal x1-x3 in die Gleichung ein um zu verifizieren, ob es eine Lösung ist.

Zitat:

da es mir nicht gelungen ist, eine Jacobimatrix numerisch zu bestimmen.

Warum willst du das denn? Die entsprechenden Löser führen so etwas intern durch.

Grüße,
Harald

Günnsen

Themenstarter

Forum-Anfänger


	Beiträge: 21

	Anmeldedatum: 26.06.16

	Wohnort: ---

	Version: ---

Verfasst am: 05.08.2016, 09:38 Titel:

Es existiert tatsächlich eine Lösung, allerdings scheint die nicht sonderlich befriedigend zu sein.

Code:

jacobi = jacobian(buy_prob);
Inverse_Matrix = inv(jacobi);
eqs_sys = prices == Inverse_Matrix * (-1 * buy_prob);

Funktion ohne Link?

Ich bin mir nicht 100% sicher, ob die Jacobi-Matrix hier bereits Teil des Lösungsverfahren ist oder nur dazu dient das GS aufzustellen.

Gruß,
Günther

Harald

Forum-Meister


	Beiträge: 24.501

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 05.08.2016, 16:08 Titel:

Hallo,

ich verstehe immer noch nicht, was dich davon abhält, fsolve ohne symbolische Variablen zu versuchen. Es ist dazu nicht nötig, selbst eine Jacobi-Matrix aufzustellen.

Grüße,
Harald

Günnsen

Themenstarter

Forum-Anfänger


	Beiträge: 21

	Anmeldedatum: 26.06.16

	Wohnort: ---

	Version: ---

Verfasst am: 05.08.2016, 18:12 Titel:

Hallo Harald,

ich befürchte, dass ich das Modell noch nicht verstanden habe. (siehe Anhang Seite 9
Formeln 4-8 )

Vielleicht siehst du ja wie ich das mit Matlab in Einklang bringen kann, weil ich sehe es nicht.

Für mich klingt das so, als ob ich Formel 7 lösen muss. Auf den Wert, der aufaddiert wird verzichte ich im ersten Schritt.

Gruß,
Günther

YAKCAY200738__akcay_et_al_paper.pdf

Beschreibung:

Download

Dateiname:

YAKCAY200738__akcay_et_al_paper.pdf

Dateigröße:

404.77 KB

Heruntergeladen:

574 mal

Günnsen

Themenstarter

Forum-Anfänger


	Beiträge: 21

	Anmeldedatum: 26.06.16

	Wohnort: ---

	Version: ---

Verfasst am: 05.08.2016, 19:49 Titel:

Hallo,

meinst du Formel 8 hart codiert liefert das Ergebnis?

Code:

function F = eqns( x )
lambda = 0.8;
theta = rand(1);
my = 1;
F(1) = lambda * exp(-(2*(x(1) - 100*theta))/my)/(my*(exp(-(x(2) - 50*theta)/my) + exp(-(x(3) - 25*theta)/my) + exp(-(x(1) - 100*theta)/my) + 1)^2) - exp(-(x(1) - 100*theta)/my)/(my*(exp(-(x(2) - 50*theta)/my) + exp(-(x(3) - 25*theta)/my) + exp(-(x(1) - 100*theta)/my) + 1))*(x(1)+80)+lambda*exp((theta*100-x(1))/my)/(1+exp((theta*100-x(1))/my)+exp((theta*50-x(2))/my)+exp((theta*25-x(3))/my));
F(2) = lambda * exp(-(2*(x(2) - 50*theta))/my)/(my*(exp(-(x(2) - 50*theta)/my) + exp(-(x(3) - 25*theta)/my) + exp(-(x(2) - 50*theta)/my) + 1)^2) - exp(-(x(2) - 50*theta)/my)/(my*(exp(-(x(2) - 50*theta)/my) + exp(-(x(3) - 25*theta)/my) + exp(-(x(2) - 50*theta)/my) + 1))*(x(2)+40)+lambda*exp((theta*50-x(2))/my)/(1+exp((theta*50-x(2))/my)+exp((theta*50-x(2))/my)+exp((theta*25-x(3))/my));
F(3) = lambda * exp(-(2*(x(3) - 25*theta))/my)/(my*(exp(-(x(2) - 50*theta)/my) + exp(-(x(3) - 25*theta)/my) + exp(-(x(3) - 25*theta)/my) + 1)^2) - exp(-(x(3) - 25*theta)/my)/(my*(exp(-(x(2) - 25*theta)/my) + exp(-(x(3) - 25*theta)/my) + exp(-(x(3) - 25*theta)/my) + 1))*(x(3)+20)+lambda*exp((theta*25-x(3))/my)/(1+exp((theta*25-x(3))/my)+exp((theta*50-x(2))/my)+exp((theta*25-x(3))/my));
end

Funktion ohne Link?

Code:

x=[0 0 0];
[result,fval,output]=fsolve(@eqns,x);
result
fval
eqns(x)

Funktion ohne Link?

Sry, dass der Code so unübersichtlich ist. Ich wollte es mal direkt ausprobieren.


Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:	Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen

Impressum | Nutzungsbedingungen | Datenschutz | FAQ |

RSS

Hosted by:

Copyright © 2007 - 2025 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.