Quy hoạch tuyến tính - Phương pháp đơn hình - C++ Code

Đây là Link video từ 1 -> 5 khá dễ hiểu đi từ cơ bản đến nâng cao về method simple trong linear programming

Tài liệu của thầy Nguyễn Văn Hiệu: Link1, Link2   (nên xem video trước khi đọc tài liệu)

Note: (cj là hệ số, xj là kết quả chúng ta cần tìm, A là bảng đơn hình)

- Dạng chuẩn tắc:

    f(x) =  sum(cj * xj)  -> min     (với j chạy từ 1 -> n, n là số ẩn x (cơ sở + ko cơ sở))

    sum(Aij * xj) = bi                 (với  i chạy từ 1 -> m, m là số ẩn x cơ sở)

    xj >= 0                                

- Cách đưa về dạng chuẩn tắc

1.    f(x) -> max    <=>    -f(x) -> min

2.    sum(Aij * xj) <= bi    <=>    sum(Aij * xj) + xNew1 = bi    

3.    sum(Aij * xj) >= bi    <=>    sum(Aij * xj) - xNew2 = bi    

4.    Nếu x ko có ràng buộc về dấu thì x = u - v     (u, v dương)

5.    ... ?

    

Problem: (bài toán đang ở dạng tổng quát)

3x1 + x2 + 3x3  -> max

2x1 +   x2 + 3x3 <= 2

  x1 + 2x2 + 3x3 <= 5

2x1 + 2x2 +   x3 <= 6

x1, x2, x3 >= 0

Solution:

B1: Nhận thấy bài toán ở dạng tổng quát nên ta đưa về dạng chuẩn tắc

-3x1 - x2 - 3x3 -> min

2x1 +   x2 + 3x3 + x4 = 2

  x1 + 2x2 + 3x3 + x5 = 5

2x1 + 2x2 +   x3 + x6 = 6

x1, x2, x3, x4, x5, x6 >= 0

B2: Lập bảng đơn hình (ở bước này ta sẽ đưa dữ liệu vào chương trình)

c = (-3, -1, -3, 0, 0, 0) // các hệ số

cb= (c4, c5, c6) //  (ở ví dụ này ta có x4 -> x5 -> x6)

=> cb = (0,0,0)

tính vector delta     //delta[i] = sum(aij * cbi) - cj

tính fx     // fx = sum(b*cb)

while(true){

    - neu vo nghiem thi ket thuc;

    - neu co nghiem thi return nghiem;

    - tìm cột vào

    - tìm cột ra

    - biến đổi ma trận sao cho cột vào thành cột ra (a, delta, b, fx biến đổi theo)

}

Code C++11: (code này được viết bởi Nguyễn Công Cường) 

#include<bits/stdc++.h>

using namespace std;

vector<int> HeSo(vector<int> c, vector< vector<double>> a){

    vector<int> cb;

    set<pair<int, int>> sp;             //dung de tim he so theo thu tu hang 1 -> hang m

   

   

    for(int j=0; j<a[0].size(); j++){

               

        int ii= -1;

        for(int i=0; i<a.size(); i++){

            if (a[i][j] == 1 && ii == -1){

                ii= i;  

            }

            else

                if (a[i][j]) ii= -2;    //co the them break de giam dpt

        }

        if (ii>=0) sp.insert({ii, j});      //neu la 1 cot trong ma tran don vi

    }

   

    for(auto p : sp)

        cb.push_back(c[p.second]);

    return cb;

}

vector<double> TinhDelta(vector<int> cb, vector<int> c, vector< vector<double>> a){

    //deltai = sum(aij * cbi) - cj

    vector<double> v;

   

    for(int j=0; j<a[0].size(); j++){

        double s= 0;

        for(int i=0; i<a.size(); i++)

            s+= a[i][j] * cb[i];

        v.push_back(s-c[j]);

    }

    return v;

}

double TinhFx(vector<int> cb, vector<double> b){

    double s= 0;

    int i= 0;

    for(auto val_b: b)

        s += val_b * cb[i++];

    return s;

}

bool VoNghiem(vector<double>delta, vector< vector<double>> a){

    //neu ton tai deltaj > 0 ma aij <= 0 (voi moi i=1:m) thi vn

   

    for(int j=0; j<delta.size(); j++)

    if (delta[j] > 0){

       

        bool check= true;

        for(int i=0; i<a.size(); i++)

            if (a[i][j] > 0) check= false;

       

        if (check) return true;

    }

    return false;

}

bool PATU(vector<double>delta){

    for(auto val:delta)

        if (val > 0) return false;

    return true;

}

vector<double> Nghiem(int m, vector< vector<double>> a, vector<double> b){

    vector<double> x;

   

    for(int j=0; j<a[0].size(); j++){

               

        int ii= -1;

        for(int i=0; i<a.size(); i++){

            if (a[i][j] == 1 && ii == -1){

                ii= i;  

            }

            else

                if (a[i][j]) ii= -2;    //co the them break de giam dpt

        }

        if (ii>=0) x.push_back(b[ii]);  //neu la 1 cot trong ma tran don vi

            else x.push_back(0);        

    }

   

    return x;

}

int TimCotVao(vector<double>delta){

    double valMax= delta[0];

    int index= 0;

   

    for(int i=1; i<delta.size(); i++){

        if (valMax < delta[i]){

            valMax= delta[i];

            index= i;

        }

    }

    return index;

}

int TimCotRa(int jVao, vector<vector<double>>a, vector<double> b, int &iVao){

    double valMin= (a[0][jVao]<=0 ? 100111000 : b[0]*1.0 / a[0][jVao]);

   

    for(int i=1; i<a.size(); i++){

       

        if (a[i][jVao] <= 0) continue;

        double tiSoKoAm= b[i]*1.0 / a[i][jVao];  

       

       

        if (valMin > tiSoKoAm){

           

            valMin= tiSoKoAm;

            iVao= i;

        }

    }

   

    for(int j=0; j<a[0].size(); j++)

    if (j != jVao){

               

        int ii= -1;

        for(int i=0; i<a.size(); i++){

            if (a[i][j] == 1 && ii == -1){

                ii= i;  

            }

            else

                if (a[i][j]) ii= -2;    //co the them break de giam dpt

        }

        if (ii == iVao) return j;  

    }

    return -1;      //error

}

void output(vector<int> c, vector<int> cb, vector<vector<double>> a, vector<double> delta, double fx, vector<double> b){

    cout<<"\n\n-----------------------BANG DON HINH--------------------------------\n\n";

    cout<<setw(10)<<""; for(auto val: c) cout<<setw(10)<<val; cout<<setw(10)<<" <-- [c]";

    cout<<"\n";

    int i= 0;

    for(auto v : a){

        cout<<"\n";

       

        cout<<setw(10)<<cb[i];

       

        for(auto val : v)

            cout<<setw(10)<<fixed<<setprecision(1)<<val;

       

        cout<<setw(10)<<b[i++];

       

    }

    cout<<"\n";

    cout<<setw(10)<<"[cb]";

    cout<<"\n";

    cout <<setw(10)<<"[delta"; for(auto val:delta) cout<<setw(10)<<fixed<<setprecision(1)<<val; cout<<"]"<<setw(10)<<fixed<<setprecision(1)<<fx<<"=fx\n";

    cout<<"\n\n-------------------------------------------------------------------------\n\n";

}

vector<double> process(int n, vector<int> c, int m, vector< vector<double>> a, vector<double> b){

   

    vector<int> cb= HeSo(c, a);

    vector<double> delta = TinhDelta(cb, c, a);

    double fx = TinhFx(cb, b);

   

    int cnt= 0;

    while(true){

        if (VoNghiem(delta, a)){

            cout<<"Bai toan vo nghiem!";

            exit(0);

        }

        if (PATU(delta)) return Nghiem(m, a, b);

       

        int jv= TimCotVao(delta);

        int iv= 0;

        int jr= TimCotRa(jv, a, b, iv);

       

        //bien doi ma tran a va vector b

        double soNhan= 1 / a[iv][jv];

        for(int j=0; j<a[0].size(); j++) a[iv][j] *= soNhan;

        b[iv] *= soNhan;

       

        for(int i=0; i<a.size(); i++)

        if (i != iv){

           

            soNhan = -a[i][jv] / a[iv][jv];    

            for(int j=0; j<a[i].size(); j++)

                a[i][j] += a[iv][j]*soNhan;

           

            b[i] += b[iv] * soNhan;

        }

       

        //bien doi vector delta

        soNhan= -delta[jv] / a[iv][jv];

        for(int j= 0; j<a[0].size(); j++)

            delta[j] += a[iv][j] * soNhan;

        fx += b[iv] * soNhan;

       

        output(c, cb, a, delta, fx, b);

    }

    return vector<double>();                                                                                              

}  

       

   

void input(){

    freopen("input.txt","r",stdin);

    //cout<<"nhap n (so an x (co so + ko co so)) : ";

    int n;  cin >> n;

   

    vector<int> c(n);  

    for(int i=0; i<n; i++){

        //cout<<"nhap c"<<i<<" : ";

        cin >> c[i];    

    }

   

    //cout<<"nhap m (so an x (co so)): ";

    int m;  cin >> m;

   

    vector< vector<double> > a(m, vector<double>(n));

    vector<double> b(m);

   

    for(int i=0; i<a.size(); i++){

        for(int j=0; j<a[i].size(); j++){

       

            //cout<<"nhap a"<<i<<j<<" : ";

            cin >> a[i][j];

        }

        //cout<<" sum(Aij * xj) = bi, nhap b"<<i<<" : ";

        cin >> b[i];

    }

   

    vector<double> x= process(n, c, m, a, b);

    cout<<"VAY NGHIEM TOI UU LAN LUOT LA: \n";

    for(int i=0; i<x.size(); i++)

        cout<<"x"<<i<<" = "<<x[i]<<"\n";

}

int main(){

   

    input();

   

}

  

Test:   (bỏ freeopen nếu nhập xuất từ console)

6

0 1 -3 0 2 0

3

1 1 -1 0 1 0 7

0 -4 4 1 0 0 12

0 -5 3 0 1 1 10  

(tren la bai tap 3 trong tài liệu thầy Hiệu)

Output bài tập 3:

Bài toán vô nghiệm

6

1 -1 0 -2 2 -2 

3

1 0 0 1 1 -1 2

0 1 0 1 0 1 12

0 0 1 2 4 3 9

(tren la bai tap 2)

Ouput của bài tập 2:

4

-5 -8

2

1 2 1 0 1

1 1 0 1 2 

(tren la bai tap 1)

4

3 -1 2 -2 

2

1 1 0 -2 1

0 -2 1 3 1

(tren la vi du 2 cua thay)

6

5 4 5 2 1 3

3

2 4 3 1 0 0 52

4 2 3 0 1 0 60

3 0 1 0 0 1 36

(tren la vi du 1 cua thay)

5

-8 0 1 1 0

3

-2 0 0 1 1 4

1 1 0 -1 0 2

1 0 1 2 0 1

(video 5)

6

-3 -1 -3 0 0 0 

3

2 1 1 1 0 0 2

1 2 3 0 1 0 5

2 2 1 0 0 1 6

(video 5)

Gradient Descent cho hàm 1 biến code C++ và Python

Problem: 

-Tìm x để f(x) min với f(x) = x^2 + 5*sin(x)

- Để giải bài toán này chúng ta có thể dùng phương pháp đã học ở cấp 3 đó là giải f'(x) = 0 và tìm cực tiểu tuy nhiên khá khó khăn 

- Một cách dễ hơn chính là dùng giải thuật Gradient Descent để có thể ước lượng được x.

Solution:

-Dự đoán một điểm x0

-Cập nhập x cho đến khi đạt điều kiện dừng

x = x0 - a * f'(x0)

x0 = x

-Điều kiện dừng:

+ Quá giới hạn số bước lặp

+ Giá trị tuyệt đối của Gradient nhỏ hơn một số gần bằng 0

+ Giá trị hàm của nghiệm tại 2 lần cập nhật không chênh lệch nhiều

Code C++:

#include<bits/stdc++.h>

#define X first

#define loop second

using namespace std;

double ff(double x){        //f'(x)

    return 2*x + 5*cos(x);

}

double f(double x){

    return pow(x, 2) + 5*sin(x);

}

typedef pair<double, int> ii;

ii GD(double alpha, double x0, double e=1e-3, int N= 1000){

   

    double x;

    for(int i=1; i<=N; i++){

       

        x = x0 - alpha * ff(x0);

        if (fabs(ff(x)) < e) return ii(x, i);

        x0 = x;

    }

    return ii(x0, N);

   

}

int main(){

    ii a= GD(0.1, -10);

    ii b= GD(0.1, 10);

    ii c= GD(0.5, 10);

    ii d= GD(0.01, 10);

   

    cout<<"Solution a.x = "<<a.X<<", f(a.X) = "<<f(a.X)<<", obtained after "<<a.loop<<" iterations\n";

    cout<<"Solution b.x = "<<b.X<<", f(b.X) = "<<f(b.X)<<", obtained after "<<b.loop<<" iterations\n";

    cout<<"Solution c.x = "<<c.X<<", f(c.X) = "<<f(c.X)<<", obtained after "<<c.loop<<" iterations\n";

    cout<<"Solution d.x = "<<d.X<<", f(d.X) = "<<f(d.X)<<", obtained after "<<d.loop<<" iterations\n";

}

Code Python:

from __future__ import division, print_function, unicode_literals

import math

import numpy as np

import matplotlib.pyplot as plt

def grad(x):

    return 2*x+ 5*np.cos(x)

def cost(x):

    return x**2 + 5*np.sin(x)

def myGD1(alpha, x0, gra = 1e-3, loop = 1000):

    x = [x0]

    for it in range(loop):

        x_new = x[-1] - alpha*grad(x[-1])

        if abs(grad(x_new)) < gra:

            break

        x.append(x_new)

    return (x, it)

if __name__ == '__main__':

    (x1, it1) = myGD1(.1, -10)

    (x2, it2) = myGD1(.1, 10)

    

    print('Solution x1 = %f, cost = %f, obtained after %d iterations'%(x1[-1], cost(x1[-1]), it1))

    print('Solution x2 = %f, cost = %f, obtained after %d iterations'%(x2[-1], cost(x2[-1]), it2))

    (x3, it3) = myGD1(.5, 10)

    print('Solution x3 = %f, cost = %f, obtained after %d iterations'%(x3[-1], cost(x3[-1]), it3))

    (x4, it4) = myGD1(.01, 10)

    print('Solution x4 = %f, cost = %f, obtained after %d iterations'%(x4[-1], cost(x4[-1]), it4))

    plt.plot(x1)

    plt.plot(x2)

    #plt.plot(x3)

    plt.plot(x4)

Output: 

Note: 

-Tốc độ hội tụ của Gradient Descent

+ Phụ thuộc vào điểm khởi tạo x0

+ Phụ thuộc vào chỉ số alpha